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1- 1  Map  of  IMS  station  locations  in  Group  1  study  region.  2 

1  -2  Study  participants  from  the  various  institutions  of  the  Group  1  IMS 

station  calibration  consortium.  3 

2- 1  Ray  path  sensitivities  calculated  using  the  extended  PL  algorithm,  (a) 

Station  DSH  and  events  for  which  sensitivities  were  calculated 
shown  with  straight  line  approximations  of  the  ray  paths,  (b),  (c), 
and  (d)  Sensitivities  for  50  km,  60  km  and  70  km  depth  slices 
calculated  for  station  DSH  and  events  shown  in  panel  (a).  Note 
that  the  sensitivities  are  unitless.  The  sensitivities  for  all  rays  that 
encounter  a  cell  are  summed,  producing  scales  that  can  range  from 
zero  upward.  Deeper  slices  through  the  sensitivity  matrix  naturally 
have  smaller  scale  ranges,  since  the  rays  spread  out  as  they 
propagate  away  from  the  station.  12 

2-2  2-D  vs.  3-D  ray  paths.  Map  view  plots  showing  the  2-D  raypath  (a) 

and  the  3-D  raypath  (b).  Depth  vs.  longitude  plots  of  the  2-D 
raypath  (c)  and  the  3-D  raypath  (d).  Travel  times  computed  by 
each  method  differ  by  0.3  sec.  In  addition  to  the  out  of  plane 
effects  seen  in  the  3-D  ray,  the  turning  depths  of  the  two  rays  differ 
by  about  5  km.  1 3 

2-3  Map  outlines  of  our  composite  3-D  velocity  model  for  the  Group  1 

study  area.  Each  of  the  submodels  has  been  defined  using  the  best 

available  reference  data  for  that  region.  The  global  background 

model  (denoted  by  the  light  blue  color)  is  the  l°-x-l°  surface  wave 

model  of  Stevens  and  Adams  (2001).  Yellow  denotes  regions 

where  the  upper  mantle  model  of  Stevens  and  Adams  has  been 

replaced  by  the  IASPEI91  upper  mantle  model  to  be  consistent 

with  local  travel  time  observations.  1 5 

2-4  Simplified  flowchart  of  the  nonlinear  joint  tomography  and  location 
procedure  used  to  develop  3-D  models  of  regions  in  Eastern  Asia 
for  the  calibration  of  Group  1  IMS  stations.  18 

2-5  While  we  only  invert  for  the  Pn  velocity,  the  velocity  profile  is 

smoothly  updated  from  the  Moho  to  410  km  based  on  the  updated 
Pn  velocity.  20 
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2-6  Results  from  a  checkerboard  resolution  test  of  the  earthquake  data  set 
used  to  update  the  WINPAK3D  model,  a)  Ray  coverage  (straight 
line  approximations)  of  the  earthquake  data  set  used  to  produce  a 
tomographic  update  to  the  Pn  velocity  over  the  model  region.  The 
blue  box  indicates  the  region  that  was  inverted  for  revised  upper 
mantle  velocities,  b)  The  model  we  used  to  predict  synthetic  travel 
times  from  events  to  stations  in  our  data  set.  c)  The  checkerboard 
pattern  resolved  by  one  iteration  of  the  nonlinear  inversion  method. 
Across  the  majority  of  the  inverted  region,  the  resolution  of  the 
individual  checkers  is  excellent.  Less  well-resolved  regions  to  the 
south-southeast  reflect  the  reduced  data  coverage,  while  poorly 
resolved  areas  in  the  far  northwest  and  southwest  reflect  sheer  lack 


of  data.  22 

2-7  Initial  (a)  and  final  (b)  Pn  velocity  maps  from  the  WINPAK3D  model 

resulting  from  application  of  our  tomography  scheme.  23 

2-8  Ray  coverage  (straight  line  approximations)  of  the  earthquakes  data  set 
used  to  produce  a  tomographic  update  to  the  Pn  velocity  over  the 
China  model  region.  25 

2-9  Initial  (a)  and  final  (b)  Pn  velocity  maps  for  the  China  region.  The  final 

Pn  velocity  model  results  from  our  nonlinear  tomography  scheme.  25 


2- 1 0  Ray  coverage  for  the  Borovoye  DSS  tomography  inversion  provided  by 

Soviet  PNE  events  (in  black)  and  Soviet  bulletin  earthquakes  (in 
red)  recorded  at  stations  of  the  Soviet  permanent  seismic  network.  26 

2- 1 1  Map  of  the  Pn  lid  velocity  corresponding  to  the  tomographically  refined 

areas  within  the  Group  1  study  region  (viz.  DSS  area  of  the  FSU, 
WINPAK  area  of  India/Pakistan,  and  the  China/Mongolia  area).  27 

3- 1  Comparison  of  P-wave  travel  time  residuals  for  55  GTO  PNEs  from  the 

FSU  measured  at  the  Borovoye  station  relative  to  IASPEI91  travel 

times  (left)  and  relative  to  the  travel  times  predicted  by  the 

regional  3-D  model  (right).  28 

3-2  Comparison  of  the  explosion  P-wave  travel  time  residuals  as  a  function 
of  distance  for  PNEs  recorded  at  the  Borovoye  station  for 
LASPEI91  travel  times  (open  circles)  and  for  the  travel  times  from 
the  3-D  regional  model  (closed  circles).  29 

3-3  Locations  of  selected  IMS  or  surrogate  stations  where  FSU  explosion 

data  have  been  used  to  validate  the  P-wave  travel  time  predictions.  29 
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3 -4a 

Comparison  of  observed  PNE  residuals  computed  with  respect  to 
IASPEI91  (left)  and  our  regional  3-D  velocity  model  (right)  for 
station  NRI. 

30 
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3 -4b 

Comparison  of  explosion  P-wave  residuals  as  a  function  of  distance  at 
NRI  computed  with  respect  to  the  IASPEI91  and  regional  3-D 
velocity  model. 
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Comparison  of  observed  PNE  residuals  computed  with  respect  to 
IASPEI91  (left)  and  our  regional  3-D  velocity  model  (right)  for 
station  TIK. 
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Comparison  of  explosion  P-wave  residuals  as  a  function  of  distance  at 
TIK  computed  with  respect  to  the  IASPEI91  and  regional  3-D 
velocity  model. 
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Comparison  of  observed  PNE  residuals  computed  with  respect  to 
IASPEI91  (left)  and  our  regional  3-D  velocity  model  (right)  for 
station  YAK. 
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Comparison  of  explosion  P-wave  residuals  as  a  function  of  distance  at 
YAK  computed  with  respect  to  the  IASPEI91  and  regional  3-D 
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Comparison  of  observed  PNE  residuals  computed  with  respect  to 
IASPEI91  (left)  and  our  regional  3-D  velocity  model  (right)  for 
station  BOD. 
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Comparison  of  explosion  P-wave  residuals  as  a  function  of  distance  at 
BOD  computed  with  respect  to  the  IASPEI91  and  regional  3-D 
velocity  model. 
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Comparison  of  observed  PNE  residuals  computed  with  respect  to 
IASPEI91  (left)  and  our  regional  3-D  velocity  model  (right)  for 
station  IRK. 
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Comparison  of  explosion  P-wave  residuals  as  a  function  of  distance  at 
IRK  computed  with  respect  to  the  LASPEI91  and  regional  3-D 
velocity  model. 
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3-9a 

Comparison  of  observed  PNE  residuals  computed  with  respect  to 
IASPEI91  (left)  and  our  regional  3-D  velocity  model  (right)  for 
station  ELT. 
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3-9b  Comparison  of  explosion  P-wave  residuals  as  a  function  of  distance  at 


ELT  computed  with  respect  to  the  IASPEI91  and  regional  3-D 
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SECTION  1 


INTRODUCTION 

1.1  BACKGROUND 

The  determination  of  accurate  seismic  locations  for  detected  events  is  one  of  the  most 
important  aspects  of  nuclear  test  monitoring  because  location  plays  such  a  key  role  in 
event  identification.  More  specifically,  a  nominal  event  location  uncertainty  area  of  1000 
km2  is  often  quoted  as  an  ultimate  location  performance  goal,  consistent,  for  example, 
with  the  onsite  inspection  area  limit  established  during  the  negotiations  leading  up  to  the 
CTBT.  However,  it  has  proven  difficult  to  demonstrate  that  this  level  of  accuracy  can  be 
routinely  achieved,  particularly  for  small  events  recorded  by  limited  numbers  of 
monitoring  stations.  The  principal  reason  for  this  difficulty  is  that  the  earth  is  laterally 
heterogeneous  over  a  very  broad  range  of  length  scales;  and,  therefore,  travel  time  curves 
based  on  radially  symmetric  earth  models,  such  as  the  standard  IASPEI91  model,  can  be 
seriously  in  error  for  some  source  and  station  locations.  This  problem  is  particularly 
acute  at  regional  distances  where  the  propagation  paths  of  the  seismic  phases  used  in 
location  are  predominantly  confined  to  the  crust  and  upper  mantle  portions  of  the  earth 
which  exhibit  the  greatest  regional  variability.  Since  many  of  the  small  seismic  events, 
which  are  of  primary  concern  in  nuclear  test  monitoring,  will  have  to  be  located  using 
regional  seismic  data,  such  variability  places  significant  limitations  on  currently 
achievable  location  accuracy.  It  follows  that  in  order  to  establish  a  robust  nuclear  test 
monitoring  capability,  it  will  be  necessary  to  calibrate  the  seismic  monitoring  stations  to 
account  for  such  systematic  deviations  from  the  nominal  travel  time  curves. 

This  report  presents  a  summary  of  the  research  which  has  been  accomplished  over  the 
past  three  years  in  a  program  directed  at  the  seismic  calibration  of  the  30  Group  1  IMS 
stations  of  Eastern  Asia  shown  in  Figure  1-1.  These  Group  1  stations  include  1 1  Primary 
and  19  Auxiliary  stations  of  the  IMS  network,  and  it  can  be  seen  from  this  figure  that  the 
area  sampled  by  the  regional  distance  ranges  around  these  stations  (i.e.  out  to  20°) 
encompasses  much  of  the  vast  landmass  of  central  and  eastern  Asia.  Simply  stated, 
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Figure  1-1.  Map  of  IMS  station  locations  in  Group  1  study  region. 

the  objective  of  our  research  has  been  to  calibrate  the  travel  time  characteristics  of  the 
crust  and  upper  mantle  structure  beneath  this  region  over  the  depth  range  sampled  by  the 
regional  seismic  ray  paths  to  these  stations  from  all  seismic  source  locations  of  potential 
interest. 

A  Consortium  of  institutions  led  by  SAIC  was  assembled  to  address  this  Group  1  IMS 
station  calibration  effort.  This  consortium  included  the  Massachusetts  Institute  of 
Technology’s  Earth  Resources  Laboratory  (ERL),  Weston  Geophysical  Corporation,  and 
the  Russian  Institute  for  Dynamics  of  the  Geospheres  (IDG).  Jack  Murphy  of  SAIC  was 
the  principal  investigator  on  the  project  and  was  responsible  for  its  overall  direction  and 
coordination.  The  MIT/ERL  group  leader  was  Nafi  Toksoz,  Jim  Lewkowicz  was  leader 
of  the  Weston  Geophysical  Corporation  effort,  and  Jamil  Sultanov  was  the  leader  of  the 
Russian  IDG  effort.  In  general,  responsibilities  for  the  different  aspects  of  this  research 
have  been  divided  between  the  participating  institutions  to  facilitate  and  focus  efforts, 
although  collaboration  between  the  groups  has  also  assisted  in  identifying  and  resolving 
specific  research  issues.  Participants  in  the  project  are  listed  in  Figure  1-2.  W.  Rodi  of 
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Figure  1-2.  Study  participants  from  the  various  institutions  of  the  Group  1  IMS  station 
calibration  consortium. 

ERL  had  overall  responsibility  for  modeling  and  statistics.  N.  Toksoz  and  his  colleagues 
at  MIT  focused  on  the  development  of  an  improved  velocity  model  for  China.  The 
Weston  group,  led  by  M.  Johnson,  had  lead  responsibility  for  ray  tracing,  implementation 
of  the  tomography  algorithm  and  computation  of  our  travel  time  correction  estimates. 

The  IDG  group,  led  by  J.  Sultanov  and  V.  Ovtchinnikov,  had  primary  responsibility  for 
preliminary  model  development  and  ground  truth  compilation  for  the  region  bounded  by 
the  territories  of  the  former  Soviet  Union. 

1.2  GENERAL  APPROACH 

Our  strategy  for  addressing  this  research  problem  was  based  on  the  belief  that,  in  order  to 
be  successful,  the  proposed  program  to  improve  regional  seismic  location  capability  in 
this  study  area  would  have  to  be  based  on  a  carefully  crafted  scientific  approach,  which 
effectively  merges  the  limited  observed  calibration  data  with  rigorous  seismological  and 
statistical  modeling  procedures.  The  ultimate  goal  is  to  generate  accurate  3-D  travel  time 
tables  for  each  station  and  seismic  phase.  Each  table  contains  predicted  travel  times  to 
the  station  from  every  point  of  a  3-D  grid  of  prospective  event  hypocenters  around  the 
station,  nominally  spaced  at  1°  distance  intervals  and  representative  depth  intervals. 
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Given  such  a  table,  the  predicted  travel  time  for  an  arbitrary  hypocenter  can  be  calculated 
via  interpolation  between  grid  points.  Subtracting  the  predicted  IASPEI91  travel  time  for 
each  grid  point  from  the  corresponding  table  travel  time  yields  a  correction  (with  respect 
to  IASPEI91)  as  a  function  of  source  hypocenter  location  for  each  IMS  station  under 
investigation.  For  the  purposes  of  this  project,  these  corrections  are  denoted  as  source 
specific  station  corrections  or  SSSCs. 

The  predicted  travel  times  for  each  station  are  represented  as  a  sum  of  model-based  travel 
times  and  empirical  corrections.  The  model-based  component  is  obtained  by  raytracing 
through  an  appropriate  3-D  velocity  model  of  the  crust  and  upper  mantle  region 
surrounding  the  station.  The  algorithm  which  we  selected  to  compute  travel  times  in  our 
complex  3-D  velocity  models  is  an  extension  of  a  finite  difference  approximation  method 
originally  developed  by  Podvin  and  Lecomte  (PL,  1991).  This  algorithm  is  well  suited  to 
the  present  application  to  rapidly  varying  models  of  the  crust  and  upper  mantle  in  that  it 
can  properly  treat  velocity  contrasts  as  great  as  10:1  independent  of  the  geometry  of  the 
feature.  The  empirical  correction  component  of  the  predicted  travel  times  is  derived  from 
calibration  events  observed  at  the  IMS  stations.  Calibration  events  include  ground  truth 
(GT)  events  with  known  locations  and  origin  times  (GTO  events),  as  well  as  GT  events 
with  various  degrees  of  uncertainty  in  their  hypocentral  coordinates  (i.e.  GT2,  GT5,  etc., 
where  the  GT  level  refers  to  the  location  uncertainty  in  km).  Thus,  the  empirical 
correction  is  essentially  an  interpolation  of  the  observed  travel  time  residuals  (with 
respect  to  our  3-D  velocity  model)  for  the  calibration  events,  analogous  to  the  empirical 
correction  surfaces  estimated  by  Schultz  et  al.  (1998)  and  others  for  different  sets  of 
stations  using  Bayesian  kriging  and  other  comparable  stochastic  inversion  procedures. 

A  key  element  of  our  approach  has  been  the  refinement  of  our  initial  3-D  velocity  model 
through  tomographic  inversion  analyses  of  relevant  data  sets,  which  have  included  the 
IMS  station  calibration  data  referenced  above  as  well  as  supplemental  data  from  the 
region  under  investigation.  These  supplemental  data  include  arrival  times  from 
calibration  events  observed  at  non-IMS  stations,  arrival  times  from  well-recorded 
earthquakes,  and  long-line  refraction  data,  where  available.  The  tomographic  inversions 
have  been  carried  out  using  a  new  algorithm  developed  for  this  project  which  is  fully 
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nonlinear  with  respect  to  both  the  velocity  model  and  event  locations,  as  is  required  for 
obtaining  accurate  solutions  in  the  structurally  heterogeneous  regional  distance  range. 
This  algorithm  has  been  implemented  as  an  iterative  nonlinear  process,  which  combines 
grid  search  event  location  and  conjugate  gradient  inversion  methods.  Within  each 
iteration,  travel  times  through  the  current  3-D  model  and  associated  sensitivities  of  travel 
times  to  cell  slownesses  are  computed  using  our  newly  extended  version  of  the  PL 
algorithm.  At  each  step  in  the  iterative  process  the  event  locations  are  re-estimated  for 
the  current  model  using  a  grid  search  location  technique  developed  by  Rodi  (2002), 
which  explicitly  accounts  for  the  varying  levels  of  GT  data  available  for  analysis. 
Following  each  event  relocation  step,  the  velocity  model  is  updated  using  a  linear 
conjugate  gradients  (LSQR)  subject  to  appropriate  smoothness  constraints.  This  entire 
process  is  then  iterated  until  the  results  converge  to  a  final  stable  estimate  of  the  3-D 
velocity  model  which  is  most  consistent  with  the  available  calibration  data.  It  is  this 
tomographically  refined  3-D  velocity  model,  together  with  associated  empirical 
corrections,  which  has  been  used  to  predict  our  final  estimates  of  the  SSSCs  for  each  of 
the  Group  1  IMS  stations,  together  with  associated  estimates  of  the  uncertainties  in  those 
corrections. 

1.3  SUMMARY  OF  MAIN  RESULTS 

As  was  noted  above,  this  report  documents  the  completion  of  a  three  year  project.  The 
main  results  from  the  studies  conducted  during  the  first  two  years  of  the  effort  have 
already  been  extensively  documented  in  previous  annual  technical  reports  (Murphy,  et 
al.,  2001,  2002)  and  will  only  be  briefly  summarized  here.  Principal  accomplishments 
achieved  under  this  project  are  described  in  more  detail  in  subsequent  sections  of  this 
report,  but  may  be  briefly  outlined  as  follows: 

•  A  new  3-D  velocity  model  of  the  entire  Group  1  study  region  has  been  formulated 
which  incorporates  both  P-  and  S-wave  velocity  distributions  defined  at  resolutions 
ranging  from  1/3  to  1  degree  across  the  area  of  interest. 

•  A  unique  Soviet  explosion  ground  truth  database  has  been  assembled  and  delivered 
to  the  SMR  (formerly  CMR)  Research  and  Development  Support  Services  (RDSS) 
data  archive.  This  database  consists  of  over  1 000  carefully  validated,  regional 
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distance  initial  P-wave  arrival  time  observations  from  Soviet  PNE  and  weapons 
tests  recorded  at  stations  of  the  FSU  permanent  seismic  network. 

•  A  sophisticated,  fully  nonlinear  tomographic  inversion  module  has  been 
formulated,  implemented,  and  applied  to  the  refinement  of  the  upper  mantle  P- 
wave  velocity  distributions  in  the  FSU,  India/Pakistan,  and  China/Mongolia 
regions  of  our  study  area.  Data  used  in  these  tomographic  inversions  have  been 
collected  into  a  uniform  database  and  delivered  to  the  RDSS  at  the  SMR. 

•  Extensive  testing  of  these  tomographically  refined  P-wave  velocity  models  against 
available  data  from  ground  truth  events  indicates  that  they  are  highly  accurate  with 
essentially  zero  bias  and  total  RMS  model  errors  on  the  order  of  1  second  across 
much  of  the  Group  1  study  area.  This  is  to  be  compared  with  an  average  bias  of 

-  3.65  seconds  and  RMS  error  of  more  than  2  seconds  associated  with  the  default 
IASPEI91  model  in  the  FSU  region. 

•  Comparison  of  regional  seismic  event  location  accuracies  for  a  sample  of  21 
widely  distributed  Soviet  PNE  explosions  using  observations  from  4-6  IMS  or 
surrogate  IMS  stations  indicate  an  average  mislocation  of  6.9  km  when  using  the 
tomographically  refined  3-D  velocity  model,  as  compared  to  an  average 
mislocation  of  20.4  km  obtained  for  the  same  data  set  using  the  default  IASPEI91 
model.  A  similar  analysis  of  data  recorded  at  a  subset  of  these  same  stations  from 
a  sample  of  20  Novaya  Zemlya  explosions  indicate  an  average  mislocation  of  less 
than  10  km  for  the  3-D  model,  as  opposed  to  an  average  mislocation  value  of  44 
km  obtained  using  the  IASPEI91  model. 

®  Similar  validation  tests  of  the  3-D  velocity  model  using  data  from  ground  truth 
events  in  the  India/Pakistan  and  China/Mongolia  regions  indicate  that  it  is  also 
significantly  superior  to  the  default  IASPEI91  model  in  those  regions.  A  ground 
truth  sample  of  arrival  time  data  observed  at  stations  of  the  Chinese  National 
Network  from  Lop  Nor  explosions  with  well-constrained  locations  (i.e.  GT1,  GT2) 
and  origin  times  has  been  collected  to  support  these  validation  studies. 

•  A  corresponding  S-wave  velocity  model  for  the  entire  Group  1  study  region  has 
been  derived  from  our  final,  tomographically  refined  P-wave  model  using 
Poisson’s  ratio  relations  determined  from  available  P  and  S  wave  velocity  data  for 
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the  crust  and  upper  mantle  across  the  area.  Limited  tests  of  this  model  have 
indicated  that  it  is  generally  consistent  with  observed  S-wave  travel  times  in  the 
region.  However,  it  is  now  clear  that  the  variance  of  the  corrected  S-wave  travel 
times  will  be  much  larger  than  that  for  the  corresponding  P-wave  times,  even  for  a 
very  good  S-wave  velocity  model,  which  means  that  they  will  be  significantly 
down-weighted  in  the  location  process. 

•  Empirical  correction  surfaces  which  account  for  travel  time  deviations  that  are  not 
accounted  for  by  the  3-D  velocity  model  have  been  estimated  using  a  newly 
developed,  multistation  kriging  algorithm.  This  new  algorithm  represents  an 
enhancement  of  the  standard  kriging  algorithms  in  that  it  explicitly  satisfies  the 
physical  constraints  of  the  seismic  reciprocity  theorem.  Application  of  these 
corrections  to  the  model-based  Soviet  explosion  ground  truth  travel  time  residuals 
results  in  reduced  RMS  errors  on  the  order  of  0.7  seconds.  We  consider  this 
residual  error  to  be  close  to  the  lower  bound  on  the  resolution  of  the  total 
calibration  process  which  we  have  implemented  for  this  study. 

•  Final  P  and  S  wave  SSSCs  have  been  estimated  for  all  30  Group  1  IMS  stations,  as 
well  as  for  1 1  surrogate  IMS  station  locations,  out  to  20  degrees  from  each  station 
for  1 1  depths  extending  from  the  earth’s  surface  down  to  200  km.  These  SSSC 
estimates  have  been  delivered  to  the  RDSS  at  SMR. 

1.4  REPORT  ORGANIZATION 

The  report  is  divided  into  seven  sections  including  these  introductory  remarks.  Section  2 
describes  the  methodology  and  the  composition  and  refinement  of  the  3-D  velocity  model 
which  form  the  bases  for  our  travel  time  prediction  methodology  and  calibration  of  the 
Group  1  study  area.  In  Section  3  we  describe  the  project  efforts  to  validate  the  travel 
time  predictions  and  to  test  refined  velocity  models  for  the  study  area.  In  Section  4  we 
discuss  extension  of  this  model-based  approach  to  travel  time  calibration  for  secondary 
seismic  phases  in  the  Group  1  study  area  and  additional  issues  arising  from  pick 
uncertainty  for  such  signals.  Section  5  describes  examples  of  the  SSSCs  derived  for  the 
regional  P-  and  S-wave  travel  times  which  have  been  derived  for  each  of  the  Group  1 
IMS  stations  and  discusses  the  significance  of  their  depth  dependence.  In  Section  6  we 
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describe  the  empirical  correction  surfaces  determined  using  the  new  reciprocal  kriging 
methodology  and  the  implications  for  location  uncertainty.  Finally,  Section  7 
summarizes  our  major  results  and  conclusions  and  addresses  issues  for  future  seismic 
travel  time  calibration  to  improve  event  location.  In  addition,  we  include  in  an  appendix 
plots  of  the  surface-focus  SSSCs  for  the  regional  P-wave  travel  times  for  each  of  the  30 
IMS  stations  from  the  Group  1  study  area. 
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SECTION  2 


METHODOLOGY  AND  VELOCITY  MODELS 
2.1  TRAVEL  TIME  PREDICTIONS 

Throughout  this  project  we  have  utilized  a  model-based  approach  to  predict  travel  times 
and  determine  SSSCs  for  the  Group  1  IMS  stations.  We  selected  this  approach  because  it 
provides  the  best  utilization  of  available  knowledge  for  use  in  regions  where  calibration 
data  are  sparse  or  non-existent.  The  essential  elements  of  this  approach  are  (1)  velocity 
models  extending  out  to  regional  distances  (i.e.  out  to  20°)  from  each  of  the  IMS  seismic 
stations  of  interest,  and  (2)  methods  of  tracing  ray  paths  and  computing  travel  times 
through  the  models  for  the  seismic  phases  useful  for  event  location.  For  the  former  we 
have  used  prior  knowledge  to  construct  a  three-dimensional  velocity  model  covering  the 
entire  Group  1  study  region,  and  we  have  used  available  ground  truth  information  to 
refine  the  model  through  an  inversion  process.  Seismic  ray  paths  through  the  model  and 
the  associated  travel  times  are  calculated  using  a  sophisticated  ray-tracing  algorithm, 
which  is  capable  of  handling  transmission  through  the  complex  3-D  structures.  The  latter 
travel  time  calculation  methods  are  used  initially  for  determining  the  residuals  forming 
the  basis  for  model  refinements  within  the  tomographic  inversion  process  and  also  for 
producing  the  travel  times  at  each  Group  1  IMS  station  from  the  regional  hypocenters 
calculated  from  the  final  3-D  velocity  model  for  eastern  Asia,  which  ultimately  determine 
the  SSSCs. 

Our  step-by-step  approach  to  the  calibration  problem  can  be  summarized  as  follows: 

(1)  Formulate  an  initial  3-D  velocity  model  for  the  Group  1  station  study  area. 

(2)  Ray  trace  through  the  initial  model  to  determine  travel  times  and  define  preliminary 
SSSCs  with  respect  to  the  IDC  default  IASPEI91  model  for  each  of  the  Group  1 
stations. 

(3)  Determine  source-specific  empirical  travel  time  corrections  for  each  station  and 
phase  by  interpolating  between  observed  calibration  event  travel  time  residuals 
relative  to  the  3-D  model. 

(4)  Generate  a  3-D  travel-time  table  for  each  Group  1  IMS  station  and  phase  by 
summing: 
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-  Model  based  travel  times  from  3-D  ray  tracing  through  the  refined  model  and 

-  Empirical  travel  time  corrections. 

(5)  Compute  revised  SSSCs  by  subtracting  IASPEI91  travel  times  from  the  resulting  3-D 
travel  time  tables  as  a  function  of  the  hypothetical  source  hypocenters  surrounding 
each  station. 

(6)  Perform  joint  tomography  and  multiple  event  location  to  update  and  refine  the 
velocity  model  and  relocate  events. 

(7)  Perform  validation  testing  of  the  refined  3-D  velocity  models  and  station  corrections 
through  relocation  analyses  of  ground  truth  events. 

(8)  Iterate  process  to  incorporate  additional  calibration  events  or  model  changes,  as  they 
become  available. 

The  method  used  to  compute  travel  times  in  the  complex  3-D  velocity  models  is  based  on 
a  finite  difference  approximation  to  the  eikonal  equation  developed  by  Podvin  and 
Lecomte  (1991)  and  implemented  by  Lomax  (1999).  To  determine  the  values  of  the 
travel  time  functions,  Tj,  we  evaluate  Tj(x;  m)  for  a  fixed  hypocenter  by  interpolating  a 
travel  time  table  stored  on  a  3-D  hypocenter  grid.  Tj  is  a  function  that  predicts  the  travel 
time  to  station  i  from  an  event  hypocenter  Xj.  This  function  depends  on  the  model 
parameter  vector  m.  The  3-D  hypocenter  grid  is  created  by  applying  the  Podvin  and 
Lecomte  (PL)  algorithm  to  the  earth  model  defined  by  m,  using  the  ith  station  as  the 
“source”  in  the  calculation.  The  PL  algorithm  can  accurately  model  different  propagation 
modes,  such  as  transmitted  and  diffracted  body  waves  or  head  waves.  The  algorithm  is 
well  suited  to  the  present  application  in  that  it  can  properly  treat  velocity  contrasts  as 
great  as  10:1,  independent  of  the  geometry  of  the  feature.  Thus,  this  approach  represents 
a  significant  improvement  over  similar  method  (e.g.  Vidale,  1988,  1990;  Moser,  1991), 
which  can  encounter  serious  difficulties  in  the  presence  of  sharp  first-order  contrasts. 

For  use  in  the  PL  algorithm,  the  model  is  discretized  on  an  equally  spaced  grid  comprised 
of  constant  velocity  cells.  The  PL  algorithm  is  used  to  compute  travel  times  from  each 
station  to  every  surrounding  grid  point.  That  is,  the  “source”  is  placed  at  the  station 
location  and,  using  reciprocity,  each  node  of  the  3-D  grid  is  then  the  hypocenter  of  an 
event.  The  resulting  grids,  one  for  each  phase  and  station,  embody  a  set  of  3-D  travel 
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time  tables  that  can  be  used  for  calculating  the  corresponding  SSSCs.  Multiple  arrivals 
(transmitted,  diffracted,  and  head  waves)  are  calculated  at  each  grid  node  and  the  first 
arrival  time  is  identified.  The  time  t  at  the  current  node  is  a  function  of  the  times  tn  at 
some  (3  or  fewer)  of  the  neighboring  nodes  and  the  slowness,  s,  in  the  cell  traversed  by 
the  wavefront  to  reach  this  node.  That  is,  t  =  t  (tn,  s).  This  method  of  travel  time 
computation  produces  a  full  grid  of  travel  times  considerably  faster  than  two-point  ray 
tracing,  and  the  sources  and  receivers  can  be  located  anywhere  within  the  model.  The  PL 
computations  are  output  in  the  form  of  3-D  travel-time  tables,  one  for  each  station  and 
seismic  phase,  which  then  can  be  used  by  a  grid-search  event  location  algorithm  in  lieu  of 
global  travel  time  tables,  such  as  IASPEI91 .  The  accuracy  of  the  PL  raytracing  algorithm 
with  respect  to  the  present  application  has  been  tested  in  a  variety  of  ways;  and  it  has 
been  concluded  that,  when  run  with  grid  spacings  on  the  order  of  5  km,  it  can  be  expected 
to  provide  3-D  travel  time  estimates  with  accuracies  of  better  than  0.25  seconds  out  to  the 
2000  km  distance  range  of  interest  in  this  calibration  study. 

The  PL  algorithm  has  also  been  extended  to  compute  the  sensitivities  of  travel  times  to 
cell  velocities,  which  are  required  as  part  of  the  model  inversion  process.  To  calculate 
the  ray  path,  we  save  the  node  pattern  (“stencil”)  at  each  step  through  the  model  as  we 
perform  the  normal  PL  forward  travel  time  calculation  to  predict  arrival  times.  This 
stencil  indicates  which  of  the  neighboring  nodes  were  used  to  calculate  the  minimum 
time  at  the  current  node.  The  stencils  can  be  used  to  reconstruct  a  path  from  any  node  of 
the  grid  to  the  source.  The  ray  tracing  is  accomplished  by  identifying  all  of  the  grid 
nodes  and  the  cells  (slownesses)  that  contribute  to  the  calculation  of  the  time  at  the 
receiver.  As  the  wavefront  propagates  away  from  the  source,  more  nodes  (and  cells)  are 
involved  in  the  travel  time  calculation  at  each  step.  After  the  midpoint  of  the  ray  path, 
the  propagation  region  narrows  until  it  reaches  the  single  node  at  the  source  location. 

The  sensitivity  of  the  travel  time  to  the  slowness,  dt/ds,  is  calculated  at  each  grid  node  of 
the  “ray”  for  the  last  cell  traversed  by  the  wavefront  to  reach  that  node.  The  weight  of 
each  neighboring  node  in  the  calculation  of  the  time  at  the  current  node  (dt/dtn) 
determines  the  weight  of  the  subsequent  node-to-source  subpath  in  the  total  travel  time 
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calculation  for  the  ray.  The  sensitivities  along  each  subpath  are  then  weighted  by  this 
term. 

This  is  illustrated  in  Figure  2-1  which  shows  the  sensitivities  of  the  travel  times 
calculated  by  PL  to  cell  slownesses  for  paths  originating  at  station  DSH  in  Tajikistan. 
After  the  sensitivities  are  calculated  for  multiple  station-centered  Cartesian  grids,  they  are 
then  mapped  to  a  single  geographic  grid.  Following  this  procedure  the  Pn  sensitivities 
are  extracted  to  form  the  kernel  matrix  for  the  tomographic  inversion  used  in  refining  the 
velocity  model.  The  ray  tracing  algorithm  and  sensitivity  calculation  were  tested  for 
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Figure  2-1 .  Ray  path  sensitivities  calculated  using  the  extended  PL  algorithm,  (a)  Station 
DSH  and  events  for  which  sensitivities  were  calculated  shown  with  straight  line 
approximations  of  the  ray  paths,  (b),  (c),  and  (d)  Sensitivities  for  50  km,  60  km  and  70  km 
depth  slices  calculated  for  station  DSH  and  events  shown  in  panel  (a).  Note  that  the 
sensitivities  are  unitless.  The  sensitivities  for  all  rays  that  encounter  a  cell  are  summed, 
producing  scales  that  can  range  from  zero  upward.  Deeper  slices  through  the  sensitivity 

matrix  naturally  have  smaller  scale  ranges,  since  the  rays  spread  out  as  they  propagate 

away  from  the  station. 
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precision  by  comparing  the  travel  time  computed  as  the  sum  of  the  weighted  sensitivity- 
slowness  products  to  the  forward  PL  calculated  times.  For  rays  with  105  nodes,  the 
difference  between  travel  times  calculated  by  these  two  methods  is  on  the  order  of  10'2 
seconds,  when  the  calculation  is  done  in  single  precision.  This  error  is  accounted  for  by 
the  level  of  precision  in  the  calculation,  demonstrating  that  the  ray  tracing  technique  is 
accurately  tracing  the  minimum  time  PL  ray  path. 


The  PL  algorithm  performs  completely  3-D  travel  time  and  ray  tracing  calculations.  It  is 
important  to  make  a  distinction  between  3-D  and  pseudo  3-D  or  “2.5-D”  methods,  where 
rays  are  traced  through  2-D  slices  of  a  3-D  model.  We  illustrate  the  difference  between 
rays  calculated  using  a  folly  3-D  method  and  a  2.5-D  method  in  an  example  of  Figure  2- 
2.  The  3-D  ray  shows  out-of-plane  effects  not  accounted  for  by  the  2-D  ray  and  a 
difference  of  5  km  in  the  turning  depth.  These  differences  arise  because  the  3-D  ray 
traverses  the  fastest  path  between  the  source  and  receiver,  which  may  or  may  not  be  in¬ 
plane;  while  the  2-D  ray  is  constrained  to  in-plane  propagation.  Consequently,  the  travel 
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Figure  2-2.  2-D  vs.  3-D  ray  paths.  Map  view  plots  showing  the  2-D  raypath  (a)  and  the 
3-D  raypath  (b).  Depth  vs.  longitude  plots  of  the  2-D  raypath  (c)  and  the  3-D  raypath 
(d).  Travel  times  computed  by  each  method  differ  by  0.3  sec.  In  addition  to  the  out  of 
plane  effects  seen  in  the  3-D  ray,  the  turning  depths  of  the  two  rays  differ  by  about  5  km, 
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time  of  the  ray  calculated  by  the  3-D  method  is  0.3  seconds  faster  than  the  travel  time  for 
the  2-D  ray.  The  discrepancy  in  travel  time  calculations  between  the  3-D  and  the  2-D 
methods  could  be  much  greater  for  even  longer  ray  paths.  We  note  that  this  example  is 
for  a  ray  with  a  source-station  distance  of  400  km,  and  our  regional  calibration  studies 
include  calculations  extending  to  20°,  or  over  2000  km  in  distance,  where  significantly 
greater  differences  could  be  expected  between  the  fully  3-D  approach  and  the  simplified 
models. 

2.2  EVENT  LOCATION 

As  noted  above,  we  incorporate  calibration  data  into  our  travel  time  calculations  and 
velocity  model  refinement  in  the  form  of  ground  truth  (GT)  events.  A  fundamental 
element  in  utilization  of  such  GT  data  is  event  location.  We  locate  events  using  the  grid 
search  method  GMEL  (Grid-Search  Multiple  Event  Location)  developed  by  Rodi  and 
Toksoz  (2002).  This  method  is  a  grid  search  technique  that  allows  for  both  Gaussian  and 
non-Gaussian  error  distributions.  GMEL  also  allows  us  to  incorporate  GT  data  of 
varying  levels,  constraining  hypocenter  solutions  based  on  the  accuracy  of  event  origins 
within  the  GT  database.  Because  GMEL  offers  the  option  of  performing  multiple  event 
locations  instead  of  traditional  single  event  locations,  we  are  able  to  use  the  calibration 
data  to  simultaneously  solve  for  a  set  of  travel  time  corrections  at  each  station.  We  do 
this  in  a  final  relocation  of  events  subsequent  to  our  tomographic  inversions  in  order  to 
account  for  any  deficiencies  in  the  crustal  part  of  our  models  that  are  not  updated  in  the 
tomographic  velocity  model  refinements. 

2.3  VELOCITY  MODELS 

2.3.1  MODEL  PARAMETERIZATION  AND  STARTING  MODELS 
We  have  defined  a  Unified  Model  Parameterization  (UMP)  that  consists  of  a  velocity- 
versus-depth  profile  at  each  point  on  a  geographic  grid,  sampled  uniformly  in  latitude  and 
longitude.  These  profiles  are  divided  into  geological  units  such  as  water,  sediment  layers, 
crustal  layers,  and  mantle  layers  that  allow  us  to  preserve  velocity  discontinuities  without 
over  parameterizing  the  model.  At  each  geographic  grid-point,  the  velocity  profile  is 
given  as  velocity/depth  pairs  at  nodes  ranging  from  sea  level  to  a  depth  of  760  km.  The 
velocity  varies  within  each  layer  as  a  smooth  gradient.  Discontinuities  in  velocity  are 
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allowed  at  the  ocean  bottom,  Moho  and  the  major  mantle  discontinuities  (410  and  660 
km  depth).  This  parameterization  reduces  the  number  of  model  parameters  to  be  solved 
for  in  the  inversion. 

Models  of  differing  parameterization  and  spatial  resolutions  are  integrated,  according  to  a 
specified  hierarchy,  and  converted  to  the  UMP  by  a  software  module  denoted  as  QUILT. 
Our  overall  3-D  velocity  model  includes  seven  submodels  covering  virtually  all  the 
Group  1  region  landmass,  plus  a  background  3-D  model  based  on  surface-wave  inversion 
analysis  (cf.  Stevens  et  al.,  2001).  The  composite  is  illustrated  in  Figure  2-3.  We 
described  the  bases  for  the  various  submodels  which  make-up  the  UMP  composite  in 
prior  annual  reports  (cf.  Murphy  et  al.,  2001, 2002).  In  defining  these  velocity 
submodels,  we  have  drawn  upon  different  knowledge  bases  depending  on  the  levels  of 
information  available  from  the  different  areas  and  utilizing  calibration  data  to  update  the 
model,  in  areas  where  that  is  appropriate.  The  initial  P-wave  velocity  model  included 
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Figure  2-3.  Map  outlines  of  our  composite  3-D  velocity  model  for  the  Group  1  study 
area.  Each  of  the  submodels  has  been  defined  using  the  best  available  reference  data  for 
that  region.  The  global  background  model  (denoted  by  the  light  blue  color)  is  the  l°-x-l° 
surface  wave  model  of  Stevens  et  al.  (2001).  Yellow  denotes  regions  where  the  upper 
mantle  model  of  Stevens  et  al.  has  been  replaced  by  the  IASPEI91  upper  mantle  model  to 
be  consistent  with  local  travel  time  observations. 
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prior  knowledge  based  on  (1)  Deep  Seismic  Sounding  results  analyzed  by  the  IDG  from 
the  FSU,  (2)  a  new  China  model  developed  by  Toksoz  and  colleagues  at  MIT  using  a 
Monte  Carlo  inversion  of  earthquake  travel  times,  (3)  a  model  for  the  India-Pakistan 
region  derived  from  published  models  and  subsequently  refined  through  tomography,  (4) 
a  new  southeast  Asia  model  based  on  previously  published  models,  (5)  models  for 
Thailand,  the  Korean  peninsula,  and  far  eastern  Russia  based  on  surface  wave  studies 
modified  by  regional  P  wave  travel  time  observations,  and  (6)  a  general  l°-x-l° 
background  model  based  on  recent  surface  wave  studies  by  Stevens  and  Adams  (2001). 
Contributions  from  the  various  models  are  assigned  priorities  and  composed  into  the 
general  3-D  velocity  model  covering  the  entire  Group  1  study  region  using  the  QUILT 
computer  software  module.  The  UMP  grids  are  generated  with  0.5°  spacing  in  latitude. 
Longitude  spacing  is  0.75°  for  the  DSS  and  Russian  tomography  regions,  while  it  is  0.5° 
for  the  WINPAK3D  and  China  tomography  regions.  This  selection  of  longitude  grid 
spacing  accounts  for  the  decreasing  aspect  ratio  of  longitude  versus  latitude  in  regions 
which  approach  the  North  Pole.  Longitude  spacing  for  the  WINPAK3D  and  China  UMP 
grids  is  defined  to  be  equal  to  the  latitude  spacing,  since  the  aspect  ratio  of  longitude 
versus  latitude  is  more  constant  over  the  latitude  ranges  for  these  model  regions. 

One  complication  in  using  this  geographically-defined  velocity  model  is  that  our  current 
implementation  of  the  Podvin-Lecomte  algorithm  solves  the  eikonal  equation  in 
Cartesian  coordinates  for  a  flat-earth  model.  To  address  this  issue,  we  have  developed 
the  algorithms  for  accurately  mapping  from  UMP  to  Cartesian  coordinates  (software 
module  XPLMOD)  as  required  by  our  ray  tracing  algorithm  and  for  mapping  the  3-D 
Cartesian  travel  time  grids  and  related  ray  path  sensitivities  back  to  geographic  grids 
(software  module  XPLTAB)  for  use  in  the  inversion.  In  mapping  from  geographic  to 
Cartesian  coordinates,  earth-flattening  transformations  are  performed  on  depth  and 
velocity  within  the  model.  A  mapping  transformation  is  carried  out  for  each  IMS  station 
to  produce  station-centered  grids.  This  preserves  distances  and  azimuths  between  all 
possible  source  locations  and  the  station,  minimizing  mapping  errors.  We  use  bi-cubic 
interpolation  to  convert  the  UMP  models  into  uniformly  sampled  Cartesian  grids.  Both 


16 


the  depths  to  layer  interfaces  and  the  layer  velocities  are  interpolated,  preserving  the 
integrity  of  geologic  layers. 

The  sampling  of  the  Cartesian  grids  is  constrained  by  limitations  of  the  finite-difference 
Podvin-Lecomte  method.  The  Podvin-Lecomte  method  uses  a  plane  wave  approximation 
which  is  inaccurate  in  the  source  region;  therefore,  the  Cartesian  grids  must  be  generated 
on  very  finely  spaced  grids  to  minimize  travel  time  error  in  this  region.  Uniformly 
sampled  5-x-5-x-5  km  spaced  grids  result  in  a  maximum  calculation  error  of  0.5  seconds 
due  to  the  plane  wave  approximation.  Using  an  enhanced  version  of  the  Podvin-Lecomte 
travel  time  algorithm  that  implements  a  multigridding  technique,  we  are  able  to  reduce 
error  in  travel  time  computation  to  less  than  0.25  seconds  at  regional  distances  for  models 
with  grid  spacing  of  5  km. 

2.3.2  MODEL  REFINEMENTS 

A  state-of-the-art  joint  velocity  tomography  and  event  relocation  algorithm  has  been 
developed  to  generate  more  accurate  velocity  models  in  geologically  complex  regions. 
Applying  this  algorithm,  we  have  generated  a  set  of  tomographically  updated  models  for 
calibration  of  the  Group  1  IMS  stations  in  three  areas:  (a)  the  India-Pakistan  region,  (b) 
China,  Mongolia,  and  surrounding  regions,  and  (c)  two  regions  in  Russia,  one 
surrounding  station  BRVK  and  the  adjacent  area  to  the  east. 

In  this  procedure  a  data  set  of  regional  seismic  events  is  used  to  iteratively  update  the 
velocity  model  following  a  conjugate  gradients  technique  that  adjusts  the  velocity  model 
to  minimize  the  misfit  between  the  calculated  and  observed  travel  times  from  multiple 
stations  and  events,  subject  to  smoothness  constraints.  The  travel  times  and  their 
sensitivities  to  the  3-D  velocity  structure  are  computed  with  an  extension  of  the  Podvin- 
Lecomte  (1991)  method.  Earthquake  locations  are  often  not  known  precisely;  therefore, 
we  also  relocate  events  using  a  3-D  grid  search  location  method  (Rodi  and  Toksoz,  2002) 
after  each  update  of  the  3-D  velocity  model.  This  process  of  inversion  is  fully  nonlinear 
with  respect  to  both  velocity  and  event  location.  However,  when  sufficient  prior  data 
exists  to  constrain  one  set  of  parameters  or  the  other,  we  can  decouple  the  location  and 
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velocity  imaging  problems.  For  example,  we  can  constrain  the  hypocenters  of  Peaceful 
Nuclear  Explosions  to  their  known  locations. 

Figure  2-4  is  a  simple  flowchart  of  the  tomography  algorithm  used  to  develop  the  new  3- 
D  velocity  model  refinements  for  the  India-Pakistan,  China-Mongolia,  and  Russian 
regions.  There  are  three  major  components  involved  in  this  joint  tomography/location 
procedure:  (1)  3-D  ray  tracing  to  predict  first  arrival  times  using  an  enhanced  version  of 
the  Podvin-Lecomte  (PL)  method  (1991);  (2)  a  3-D  grid  search  location  algorithm 
(GMEL)  by  Rodi  and  Toksoz  (2002)  to  relocate  events  inside  the  appropriate  velocity 
model;  and  (3)  a  linear  conjugate  gradient  inversion  algorithm  to  produce  the  updated 
velocity  model  inside  each  iteration  of  the  overall  nonlinear  algorithm.  In  the  next  few 
subsections,  we  describe  the  tomographic  velocity  model  refinements  which  make-up  the 
new  3-D  velocity  model  for  the  Group  1  region. 


Figure  2-4.  Simplified  flowchart  of  the  nonlinear  joint  tomography  and  location 
procedure  used  to  develop  3-D  models  of  regions  in  Eastern  AsiaTorlhe_calibration-Qf 

Group  1  IMS  stations. 
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In  the  preceding  subsections,  we  described  the  3-D  raytracing  using  the  PL  algorithm  and 
utilization  of  the  GMEL  for  relocating  events  within  the  region  covered  by  the  model. 
Once  the  database  of  events  has  been  relocated  and  the  ray  path  sensitivities  have  been 
calculated,  we  perform  a  linear  conjugate  gradients  inversion  for  an  optimized  update  to 
the  velocity  model.  We  use  a  version  of  the  LSQR  algorithm  (cf.  Nolet,  1983;  Paige  and 
Saunders,  1982)  to  produce  this  update.  The  LSQR  algorithm  is  a  linear  conjugate 
gradient  method  used  to  iteratively  solve  large  systems  of  sparse,  linear  equations.  The 
output  of  the  algorithm  is  a  vector  of  changes  to  the  input  velocity  model.  The  model 
update  produced  by  LSQR  is  constrained  in  several  ways.  First,  we  fix  a  small  buffer 
region  along  the  model  region  perimeter  to  the  values  of  the  initial  model.  This  is  to 
prevent  large  velocity  variations  from  occurring  in  areas  of  the  model  that  are  poorly 
constrained  by  data  and  to  allow  us  to  seamlessly  integrate  our  final  models  into  other 
global  models.  Second,  we  apply  a  smoothing  operator  to  the  model  using  a  second 
differencing  operator,  which  is  equivalent  to  ensuring  the  curvature  of  the  model  is 
smooth  (Twomey,  1977).  Finally,  we  apply  a  scalar  damping  parameter  to  the  model  to 
balance  the  sharpness  or  noisiness  with  the  horizontal  spread  or  smoothness  of  the 
recovered  velocity  contrasts.  After  the  linear  conjugate  gradient  method  has  converged 
to  an  optimized  update  to  the  model,  we  use  the  model  change  vector  as  a  search 
direction  in  the  next  iteration  of  the  nonlinear  conjugate  gradient  inversion. 

Thus,  our  inversion  algorithm  jointly  solves  the  nonlinear  problem,  iterating  over  linear 
inversion  steps  that  include  updates  of  the  hypocenters,  velocities  and  ray  paths.  This 
technique  explicitly  addresses  the  coupling  between  the  hypocenters  and  the  velocity 
structure  by  computationally  breaking  down  the  large  matrix  that  must  be  inverted  for 
velocities  and  locations.  This  computational  technique  results  in  two  separate  smaller 
matrices  which  may  be  inverted  separately,  but  still  solve  the  simultaneous  problem  (cf. 
Spencer  and  Gubbins,  1980;  Rodi  et  al.,  1981). 

Although  we  are  currently  only  solving  for  Pn  velocity  in  the  inversion,  model  updates 
extend  beyond  Pn  velocity.  Based  on  the  updated  Pn  velocity,  velocities  from  the  Moho 
to  410  km  are  refined  in  a  smooth  manner  as  shown  in  Figure  2-5.  A  constant  shift  in 
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Figure  2-5.  While  we  only  invert  for  the  Pn  velocity,  the  velocity  profile  is 
smoothly  updated  from  the  Moho  to  410  km  based  on  the  updated  Pn  velocity. 

velocity  is  applied  from  the  Moho  extending  to  2 10  km  depth.  Tapering  between  210  km 
and  410  km  depth  insures  a  smooth  transition  down  to  the  4 10  km  discontinuity. 

Application  to  the  India-Pakistan  Region 

We  first  look  at  tomography  applied  in  the  India-Pakistan  region  using  a  set  of 
earthquake  data.  The  Weston  Geophysical  India-Pakistan  3-D  velocity  model 
(WINPAK3D),  a  detailed  3-D  velocity  model  synthesized  from  numerous  previous 
studies,  was  used  as  an  initial  model.  The  references  utilized  to  develop  the  initial 
WTNPAK3D  model  included  data  such  as  seismic  reflection  and  refraction  surveys  (i.e. 
DSS  profiles,  Pn  tomography,  Pnl  waveform  inversion),  interpretations  of  gravity  data, 
surface  wave  studies,  and  receiver  function  analyses.  The  velocity  model  is  defined  on  a 
grid  of  one-degree  by  one-degree  blocks  and  5  km  depth  intervals  from  0  to  75  km.  The 
IASPEI91  model  (Kennett  and  Engdahl,  1991)  is  appended  to  the  base  of  the  preliminary 
velocity  model,  beginning  at  80  km  depth  and  extending  into  the  mantle,  to  accommodate 
ray  paths  that  travel  into  the  upper  mantle.  Some  preliminary  validation  was  performed 
on  the  initial  model  to  verify  its  suitability  and  potential  to  improve  event  locations.  See 
Johnson  and  Vincent  (2002)  for  details  on  the  development  and  validation  of  this  initial 
WINPAK3D  model.  Previously,  the  eastern  boundary  of  the  India-Pakistan  model,  at 
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80°  longitude,  encompassed  only  central  and  western  India.  Recently,  Weston 
Geophysical  extended  the  bounds  of  the  WINPAK3D  region  to  incorporate  all  of  India 
and  performed  joint  velocity  and  hypocenter  relocation  tomography  on  the  new  model, 
solving  for  revised  upper  mantle  velocities.  This  work  was  completed  under  the  auspices 
of  the  Defense  Threat  Reduction  Agency  contract  number  DTRA01-00-C-0098. 

The  data  set  for  the  inversion  is  made  up  of  a  larger,  more  accurate  database  of  events 
than  was  used  in  the  initial  WINPAK3D  inversion.  Events  were  selected  from  an 
updated  version  of  the  Engdahl  et  al.  (1998)  database  (EHB)  of  well-located  earthquakes 
and  explosions  that  we  obtained  in  February  2002.  The  EHB  data  set  is  generally 
considered  to  be  GT15  or  better  in  continental  areas,  according  to  the  IASPEI  Working 
Group  on  Reference  Events  (http://lemond.colorado.edu/~copgte/').  To  ensure  that  the 
most  accurate  data  subset  was  used  in  the  inversion,  only  events  that  had  a  secondary 
station  azimuth  gap  (largest  azimuthal  gap  closed  by  one  station)  of  <130°  were  selected 
from  the  EHB  data  set.  Based  on  relocated  explosion  studies,  95%  of  the  events  will  be 
mislocated  by  no  worse  than  15  km  (Engdahl,  2002).  Additionally,  data  outliers  with 
residuals  of  greater  than  ±  7  sec  were  omitted.  This  threshold  was  chosen  because 
residuals  of  5  seconds  or  greater  can  be  attributed  to  structure  in  tectonically  complicated 
regions  such  as  the  Hindu  Kush  (cf.  Murphy  et  al.,  1998).  We  chose  a  distance  cut-off  of 
14°  in  order  to  filter  out  rays  bottoming  deeper  than  the  upper  mantle.  We  also  restricted 
event  depths  to  the  crust  and  upper  mantle,  with  a  maximum  depth  of  60  km.  Final 
selection  was  focused  on  obtaining  optimal  spatial  distribution  across  the  region,  both  in 
latitude/longitude  and  depth.  The  events  in  the  data  set  were  recorded  at  89  stations 
within  our  model  region,  and  consisted  of  over  22,500  arrivals  for  1,337  events  that  were 
used  in  the  joint  inversion.  This  data  set  is  shown  in  Figure  2-6a. 

To  ensure  that  the  tomography  algorithm  was  functioning  properly  and  to  provide  a  test 
of  the  resolution  of  our  data  set  for  the  India-Pakistan  region,  we  performed  a 
“checkerboard”  inversion  test  using  the  subset  of  the  EHB  database  (Figure  2-6a). 

For  the  resolution  test,  we  perturbed  the  Pn  velocity  of  an  IASPEI91  1-D  model, 
imposing  a  3°  by  3°  checkerboard  pattern  of  alternating  high  (viz.  8.3  km/sec)  and  low 
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Figure  2-6.  Results  from  a  checkerboard  resolution  test  of  the  earthquake  data  set  used 
to  update  the  WINPAK3D  model,  a)  Ray  coverage  (straight  line  approximations)  of  the 
earthquake  data  set  used  to  produce  a  tomographic  update  to  the  Pn  velocity  over  the 
model  region.  The  blue  box  indicates  the  region  that  was  inverted  for  revised  upper 
mantle  velocities,  b)  The  model  we  used  to  predict  synthetic  travel  times  from  events  to 
stations  in  our  data  set.  c)  The  checkerboard  pattern  resolved  by  one  iteration  of  the 
nonlinear  inversion  method.  Across  the  majority  of  the  inverted  region,  the  resolution 
of  the  individual  checkers  is  excellent.  Less  well-resolved  regions  to  the  south- 
southeast  reflect  the  reduced  data  coverage,  while  poorly  resolved  areas  in  the  far 
northwest  and  southwest  reflect  sheer  lack  of  data. 


(viz.  7.9  km/sec)  velocities  (cf.  Figure  2-6b).  Inversion  was  performed  using  synthetic 
arrival  times  calculated  through  the  perturbed  IASPEI91  model,  while  the  unmodified 
LASPEI91  model  (Pn  velocity  of  8.04  km/s)  served  as  the  starting  model.  We  performed 
one  iteration  of  our  nonlinear  conjugate  gradient  scheme  to  retrieve  an  estimate  of  the 
resolving  capability  of  our  data  set.  The  damping  parameter  was  chosen  to  reduce  the 
RMS  error  while  keeping  the  noise  (one-node  variations)  low.  The  damping  parameter 
chosen  also  preserves  the  amplitude  of  the  velocity  variations  of  the  synthetic  model. 
Figure  2-6c  shows  that  the  checkerboard  pattern  Pn  velocity  was  successfully  recovered 
in  the  inversion  with  especially  good  results  in  the  central  areas  having  very  dense  ray 
coverage.  These  results  indicate  both  that  our  tomography  algorithm  is  performing  well 


and  that  the  set  of  earthquake  data  selected  to  update  the  WINPAK3D  model  is  providing 
high.reso1ution-across-much-of-theL.reg-ion, - 
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After  the  checkerboard  resolution  test,  we  performed  several  iterations  of  our  nonlinear 
inversion  technique  to  update  the  Pn  velocity  map  extracted  from  our  initial  model.  We 
inverted  for  Pn  velocity  as  a  function  of  latitude  and  longitude  and  then  imposed  a  fixed 
rule  for  extrapolating  this  velocity  into  the  upper  mantle.  Tomography  was  applied  to  the 
initial  WINPAK3D  model  using  the  arrivals  shown  in  Figure  2-6a.  The  event 
hypocenters  were  constrained  to  the  EHB  locations,  since  the  starting  locations  were 
known  to  be  very  accurate  based  on  the  event  selection  criteria.  The  damping  parameter 
for  the  inversions  was  again  chosen  to  reduce  the  RMS  while  keeping  the  noisiness  of  the 
recovered  velocity  change  low.  Figure  2-7a  shows  the  initial  Pn  map.  The  final  Pn 
model  is  presented  in  Figure  2-7b.  The  general  distribution  of  high  and  low  velocities  in 
the  final  model  is  similar  to  the  starting  model,  but  contains  more  detail.  The  RMS 
travel-time  error  of  the  tomography  data  set  was  1 .70  sec  for  the  updated  model, 
compared  to  2.28  sec  for  the  IASPEI91  model. 


Velocity  (km/s)  Velocity  (km/s) 


Figure  2-7.  Initial  (a)  and  final  (b)  Pn  velocity  maps  from  the  WINPAK3D  model 
resulting  from  application  of  our  tomography  scheme. 
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Application  to  the  China  Region 

The  preliminary  model  for  the  China  inversion  was  developed  at  M.I.T.  (Toksoz  et  ah, 
2003)  as  a  composite  of  I  D  models  defined  on  a  latitude-longitude  grid,  with  the  model 
at  each  grid  point  being  obtained  by  Monte  Carlo  inversion  of  P-wave  travel-times  from 
rays  lying  within  a  10-degree  square  region  centered  on  its  grid  point.  Development  of 
the  initial  China  model  was  described  in  detail  in  our  prior  annual  report.  We  updated 
this  model  using  a  subset  of  earthquake  data  selected  from  the  Annual  Bulletin  of 
Chinese  Earthquakes  (ABCE)  that  were  chosen  to  ensure  a  reliable  data  set  and  optimum 
coverage  for  tomographic  inversion.  The  criteria  for  event  selection  were  similar  to  those 
used  for  the  WINPAK3D  database.  Data  outliers  with  residuals  of  greater  than  ±  7  sec 
were  omitted.  We  again  chose  a  distance  cut-off  of  14°  in  order  to  filter  out  rays 
bottoming  deeper  than  the  upper  mantle.  We  also  restricted  event  depths  to  the  crust  and 
upper  mantle,  with  a  maximum  depth  of  60  km.  The  resulting  database  included  over 
64,000  arrivals  from  nearly  10,000  events  recorded  at  88  stations  in  and  around  China 
(Figure  2-8). 

We  inverted  for  Pn  velocity  as  a  function  of  latitude  and  longitude  and  extrapolated  this 
velocity  into  the  upper  mantle  as  described  previously.  Tomography  was  applied  to  the 
initial  China  model  using  the  arrivals  shown  in  Figure  2-8.  In  this  case,  we  performed 
hypocenter  relocation  as  part  of  the  tomographic  inversion  procedure.  Events  were 
relocated  with  varying  constraints  relative  to  the  initial  locations  provided  in  the  ABCE 
bulletin.  The  new  hypocenter  solutions  were  constrained  to  be  within  a  maximum  of  1 5 
km  from  the  starting  location,  in  epicenter  and  in  depth.  The  revised  model  and  event 
locations  reduced  the  root  mean  square  (RMS)  travel-time  residuals  as  compared  with  the 
global  1-D  IASPEI91  model  from  2.84  sec  to  1.55  sec.  The  initial  velocity  model  (Figure 
2-9a)  is  less  variable  in  velocity  compared  with  the  updated  Pn  model  (Figure  2-9b).  The 
inverted  model  also  appears  to  resolve  model  details  with  higher  resolution. 

Application  to  the  Former  Soviet  Union 

We  have  described  the  tomographic  inversion  procedures  applied  to  areas  of  the  FSU 
with  good  groundtruth  data  in  prior  annual  reports.  In  particular,  tomographic  analyses 
similar  to  those  for  the  WINPAK  region  were  applied  to  the  FSU  region  around  IMS 
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Figure  2-8.  Ray  coverage  (straight  line  approximations)  of  the  earthquake  data  set 
used  to  produce  a  tomographic  update  to  the  Pn  velocity  over  the  China  model  region. 


Figure  2-9.  Initial  (a)  and  final  (b)  Pn  velocity  maps  for  the  China  region.  The  final  Pn 
velocity  model  results  from  our  nonlinear  tomography  scheme. 


25 


station  BRV  and  the  adjacent  region  to  the  east  in  Russia.  For  these  analyses  we  used 
data  from  PNE  explosions  throughout  the  FSU  plus  additional  data  from  nuclear  tests  at 
the  Semipalatinsk  and  Novaya  Zemlya  test  sites.  These  explosions  were  recorded  at 
BRV  and  other  stations  of  the  FSU  permanent  seismic  network,  including  several  with 
locations  at  or  near  IMS  station  sites.  For  the  Borovoye  station  the  explosion  data  were 
supplemented  with  about  300  precise  arrival  time  observations  at  BRVK  from  Soviet 
bulletin  earthquakes,  located  primarily  along  the  active  seismic  zones  to  the  south  and 
east.  Figure  2-10  illustrates  travel  paths  and  ray  coverage  used  in  the  tomographic 
inversion  for  the  region  around  BRVK  for  FSU  nuclear  explosions  and  the  supplemental 
earthquakes.  So,  over  much  of  this  FSU  region  we  have  very  accurate  travel  times  and 
good  station  coverage  from  the  PNE  database  which  enables  very  good  resolution  for  the 
tomographic  inversion. 


Figure  2-10.  Ray  coverage  for  the  Borovoye  DSS  tomography  inversion  provided 
by  Soviet  PNE  events  (in  black)  and  Soviet  bulletin  earthquakes  (in  red)  recorded 
at  stations  of  the  Soviet  permanent  seismic  network. 
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As  with  the  other  two  subareas,  we  found  that  the  tomographically  revised  models  for  the 
FSU  regions  produced  reductions  in  the  travel-time  errors  relative  to  IASPEI91 .  The 
RMS  travel-time  error  for  the  Soviet  PNEs  from  the  tomography  data  set  was  1.1  sec  for 
the  updated  model,  compared  to  3.1  sec  for  the  IASPEI91  model.  Results  of  the 
tomographic  analyses,  including  the  FSU,  WINPAK,  and  China  region  are  shown  in 
Figure  2-11.  We  show  here  the  Pn  velocities  for  the  tomographically  refined  models 
covering  much  of  our  Group  1  region.  Focusing  on  the  FSU  region,  we  find  some 
apparent  correlations  of  the  tomography  model  results  with  published  crust  and  upper 
mantle  characteristics  in  the  region.  For  example,  two  lower  velocity  features  to  the 
north  correlate  remarkably  well  with  regions  of  crustal  thinning  associated  with  the  NE 
Volga-Ural  Uplift  (to  the  west)  and  the  crest  of  the  West  Siberian  Platform  (to  the  east). 
Similarly,  a  prominent  low  velocity  zone  in  eastern  China  correlates  with  the  Yin  Shan 
Uplift.  So,  it  appears  that  the  tomographic  analyses  are  generally  giving  reasonable 
results. 
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Figure  2-11.  Map  of  the  Pn  lid  velocity  corresponding  to  the  tomographically 
refined  areas  within  the  Group  1  study  region  (viz.  DSS  area  of  the  FSU, 
WINPAK  area  of  India/Pakistan,  and  the  China/Mongolia  area). 
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SECTION  3 


TRAVEL  TIME  PREDICTION  VALIDATION 
3 . 1  FORMER  SOVIET  UNION  VALIDATION 

Our  validation  analyses  of  the  velocity  models  for  the  FSU  have  been  described 
extensively  in  prior  annual  reports;  we  show  in  this  section  some  of  the  highlights  of 
those  prior  reports  as  well  as  some  additional  recent  developments.  Figure  3-1  shows  a 
comparison  of  the  P-wave  travel-time  residuals  for  55  FSU  PNEs  measured  at  the 
Borovoye  station  (BRV)  relative  to  the  nominal  IASPEI91  travel  times  (left)  and  relative 
to  our  3-D  velocity  model  (right).  It  is  apparent  from  the  left-hand  figure  that  the 
observed  travel  times  are  much  faster  than  IASPEI91  for  most  of  the  events,  as  evidenced 
by  the  large  negative  residuals.  On  the  other  hand,  for  our  3-D  regional  velocity  model 
the  travel  time  residuals  are  significantly  reduced  and  are  more  random  geographically. 

Figure  3-2  shows  the  same  residuals  at  BRV  plotted  as  a  function  of  event  distance  for 
the  FSU  PNEs.  The  average  residual  relative  to  IASPEI91  travel  times  is  -  3.5  seconds 
compared  to  average  residuals  of  only  0.06  seconds  for  the  3-D  model.  Furthermore, 
there  is  also  a  significant  reduction  in  the  scatter  of  the  data  from  a  =  2.26  seconds  for 
IASPEI91  to  cr  -  0.96  seconds  for  the  3-D  model. 
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Figure  3-1.  Comparison  of  P-wave  travel  time  residuals  for  55  GT0  PNEs  from  the 
FSU  measured  at  the  Borovoye  station  relative  to  IASPEI91  travel  times  (left)  and 
relative  to  the  travel  times  predicted  by  the  regional  3-D  model  (right). 
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DISTANCE 

Figure  3-2.  Comparison  of  the  explosion  P-wave  travel  time  residuals  as  a  function 
of  distance  for  PNEs  recorded  at  the  Borovoye  station  for  IASPEI91  travel  times 
(open  circles)  and  for  the  travel  times  from  the  3-D  regional  model  (closed  circles). 

We  made  similar  comparisons  for  the  PNEs  recorded  at  8  additional  IMS  and  surrogate 
station  sites  in  the  FSU  (cf.  Figure  3-3)  with  similar  results.  For  each  of  the  stations,  we 
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Figure  3-3.  Locations  of  selected  IMS  or  surrogate  stations  where  FSU 
explosion  data  have  been  used  to  validate  the  P-wave  travel  time  predictions. 
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computed  the  observed  P-wave  travel  time  residuals  with  respect  to  the  predictions  based 
on  the  IASPEI91  model  and  our  refined  3-D  velocity  model  for  the  region.  Figures  3-4 
through  3-11  show  the  P-wave  residuals  plotted  as  a  function  of  location  and  distance 
based  on  the  FSU  PNE  ground  truth  events.  The  results  indicate  that  in  all  cases  the 
residuals  are  consistently  reduced  for  our  3-D  model  relative  to  the  IASPEI91  model. 
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Figure  3-4a.  Comparison  of  observed  PNE  residuals  computed  with  respect  to 
IASPEI91  (left)  and  our  regional  3-D  velocity  model  (right)  for  station  NRI. 
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Figure  3-4b.  Comparison  of  explosion  P-wave  residuals  as  a  function  of  distance 
at  NRI  computed  with  respect  to  the  IASPEI91  and  regional  3-D  velocity  model. 
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Figure  3-5a.  Comparison  of  observed  FNE  residuals  computed  with  respect  to 
IASPEI91  (left)  and  our  regional  3-D  velocity  model  (right)  for  station  TIK. 
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Figure  3-5b.  Comparison  of  explosion  P-wave  residuals  as  a  function  of  distance  at 
TIK  computed  with  respect  to  the  IASPEI91  and  regional  3-D  velocity  model 
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Figure  3 -6a  Comparison  of  observed  PNE  residuals  computed  with  respect  to 
IASPEI91  (left)  and  our  regional  3-D  velocity  model  (right)  for  station  YAK. 
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Figure  3-6b.  Comparison  of  explosion  P-wave  residuals  as  a  function  of  distance  at 
YAK  computed  with  respect  to  the  IASPEI91  and  regional  3-D  velocity  model. 
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Figure  3-7a.  Comparison  of  observed  PNE  residuals  computed  with  respect  to 
IASPEI91  (left)  and  our  regional  3-D  velocity  model  (right)  for  station  BOD. 
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Figure  3-7b.  Comparison  of  explosion  P-wave  residuals  as  a  function  of  distance 
at  BOD  computed  with  respect  to  the  IASPEI91  and  regional  3-D  velocity  model. 
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Figure  3-8a.  Comparison  of  observed  PNE  residuals  computed  with  respect  to 
IASPEI91  (left)  and  our  regional  3-D  velocity  model  (right)  for  station  IRK.. 
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Figure  3-8b.  Comparison  of  explosion  P-wave  residuals  as  a  function  of  distance  at 
IRK  computed  with  respect  to  the  IASPEI91  and  regional  3-D  velocity  model. 
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Figure  3-9a.  Comparison  of  observed  PNE  residuals  computed  with  respect  to 
I  ASPEN  1  (left)  and  our  regional  3-D  velocity  model  (right)  for  station  ELT. 
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Figure  3-9b.  Comparison  of  explosion  P-wave  residuals  as  a  function  of  distance  at 
ELT  computed  with  respect  to  the  IASPEI91  and  regional  3-D  velocity  model. 
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Figure  3- 10a.  Comparison  of  observed  PNE  residuals  computed  with  respect  to 
IASPEI91  (left)  and  our  regional  3-D  velocity  model  (right)  for  station  FRU. 


at  FRU  computed  with  respect  to  the  IASPEI91  and  regional  3-D  velocity  model. 
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Figure  3-1  la.  Comparison  of  observed  PNE  residuals  computed  with  respect  to 
LA.SPEI91  (left)  and  our  regional  3-D  velocity  model  (right)  for  station  Ml  1 . 
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Figure  3-1  lb.  Comparison  of  explosion  P-wave  residuals  as  a  function  of  distance 
at  Ml  1  computed  with  respect  to  the  IASPEI91  and  regional  3-D  velocity  model. 


As  noted  above,  in  all  cases  the  average  travel  time  residuals  have  been  greatly  reduced 
using  the  3-D  regional  model  (cf.  Table  3-1).  Over  the  9  stations  the  average  P-wave 
residual  bias  drops  from  -  3.65  seconds  for  IASPEI91  to  -  0.10  seconds  for  the  3-D 

model.  The  corresponding  average  scatter  is  reduced  from  o  =  2.12  seconds  for 

IASPEI91  to  o  =  0.94  seconds  for  the  3-D  model.  So,  our  3-D  velocity  model  appears 
to  be  very  accurately  representing  the  travel  times  throughout  the  DSS  model  region  of 
the  FSU  sampled  by  the  widely  distributed  Soviet  PNE  events. 


Table  3-1 .  Comparison  of  average  P-wave  travel  time  residuals  for  PNEs  recorded 
at  IMS  stations  in  the  FSU  based  on  IASPEI91  and  3-D  regional  model. 

Summary  Soviet  PNE  Statistics  for 
Tomographicallv  Refined  DSS  Velocity  Model 

IASPEI-91  3-D  Model 


Station 

At,  seconds 

a, seconds 

At,  seconds 

a, seconds 

BRVK 

-3.50 

2.26 

0.06 

0.96 

NRI 

-3.34 

2.93 

-0.12 

0.98 

TIK 

-4.33 

1.38 

0.27 

0.90 

YAK 

-3.80 

2.43 

0.38 

1.23 

BOD  (PDY) 

-4.26 

2.52 

-0.56 

0.93 

IRK  (TLY) 

-3.46 

2.67 

-0.48 

0.85 

ELT  (ZAL) 

-4.10 

1.66 

-0.03 

0.74 

FRU  (AAK) 

-2.06 

1.36 

-0.07 

1.14 

Mil  (OBN) 

-4.01 

1.83 

-0.31 

0.77 

At  =  -3.65  sec 

~o-  2.12  sec 

At  =  -0.10  sec 

o’-  0.94  sec 

To  further  assess  the  relevance  of  the  improvement  to  the  P-wave  travel  time  estimates 
produced  from  our  3-D  regional  model,  we  have  analyzed  the  P-wave  errors  associated 
with  various  travel  time  representations  used  in  event  location.  In  particular,  in  Figure  3- 
12  we  compare  the  travel  time  residual  errors  as  a  function  of  distance  from  our 
tomographically  refined  3-D  model  for  the  DSS  region  with  several  alternative  P-wave 
error  estimates.  The  “Original  IDC”  error  corresponds  to  the  model  used  in  weighting 
the  travel  time  observations  within  the  IDC  hypocenter  location  code;  so,  it  essentially 
corresponds  to  the  errors  relative  to  the  default  IASPEI91  travel  time  tables  and  is  very 
conservative.  The  error  model  labeled  “CUB”  corresponds  to  results  from  the  Group  2 
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Comparison  of  P-Wave  Misfit  Statistics  for  Various  Velocity  Models  , 


Figure  3-12.  P-wave  error  estimates  as  a  function  of  distance  based  on  travel  time 
residuals  from  our  3-D  models  for  the  DSS  region  compared  to  comparable  errors 
from  other  regions. 

calibration  effort  and  is  based  on  the  residuals  for  events  in  the  Group  2  region  relative  to 
the  CUB  velocity  model.  Our  “Original  DSS”  model  is  a  summary  based  on  the  travel 
time  residuals  for  our  PNE  ground  truth  dataset  of  about  1000  P-wave  observations 
relative  to  our  initial  DSS  velocity  model;  and  the  “Tomographically  Refined  DSS”  curve 
is  a  similar  plot  for  the  residuals  after  tomographic  refinements  to  the  3-D  model.  In 
general,  the  three  regional  models  all  show  reduced  error  in  the  regional  distance  range 
relative  to  the  errors  from  the  IASPEI91  model  (i.e.  Original  IDC).  For  our  Original  DSS 
model  the  associated  travel  time  errors  show  no  particular  dependence  on  distance  over 
this  regional  distance  range,  and  the  overall  level  is  comparable  to  that  for  the 
uncalibrated  teleseismic  travel  times,  which  is  basically  the  level  of  the  Original  IDC 
curve  over  the  flat  portion  for  A  >  25°.  So,  considering  the  comparable  levels  of  the 
errors,  the  regional  P-wave  travel  times  with  the  corrections  from  the  Original  DSS 
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model  should  be  equally  weighted  with  the  teleseismic  observations  in  seismic  event 
location  algorithms.  Furthermore,  the  corresponding  statistics  for  our  Tomographically 
Refined  DSS  model  show  again  essentially  no  distance  dependence  and  an  average  c 
level  substantially  lower  than  the  nominal  teleseismic  ct  (i.e.  about  1 .0  seconds  or  less  for 
the  refined  regional  model  versus  ~1 .4  seconds  for  teleseismic).  It  should  be  noted  that 
this  result  may  be  a  little  too  optimistic,  since  the  data  from  which  the  errors  were  derived 
were  also  used  in  the  tomographic  refinement  process.  However,  it  seems  likely  that  the 
true  travel  time  prediction  error  from  our  regional  model  for  the  DSS  region  of  the  FSU 
should  lie  between  the  “Original  DSS”  and  the  “Tomographically  Refined  DSS”  curves. 
In  any  case,  there  seems  to  be  no  strong  evidence  of  distance  dependence  in  the  regional 
P-wave  travel  time  error  estimates  from  our  regional  3-D  model,  and  the  a  values  are 
comparable  to  or  less  than  those  for  the  uncalibrated  teleseismic  observations.  So,  it 
appears  that  application  of  these  regional  travel  time  corrections  should  lead  to  better 
determination  of  seismic  event  locations. 

Of  course,  a  critical  test  of  performance  for  the  3-D  velocity  model  for  events  from  the 
FSU  is  the  accuracy  of  event  location.  We  analyzed  the  location  accuracy  for  14  GTO 
PNE  events  from  the  FSU  recorded  by  4  to  6  regional  IMS  stations  (cf.  Figure  3-13). 
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Figure  3-13.  Map  locations  of  14  FSU  PNE  events  (circles)  used  in  analyzing 

location  accuracy  for  the  3-D  regional  velocity  model.  Each  explosion  was 
recorded  at  4-6  regional  IMS  station  locations  (triangles)  in  the  Group  1  study  area. 
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Figure  3-14  shows  a  comparison  of  the  location  errors  for  these  events  based  on  locations 
using  nominal  IASPEI91  travel  times  versus  the  locations  using  our  3-D  velocity  model 
for  the  region.  We  see  that  in  nearly  all  cases  the  mislocations  are  reduced  using  the  3-D 
model.  For  the  14  PNEs  the  average  mislocation  for  the  3-D  model  is  only  6.9  km 
compared  to  an  average  mislocation  of  20.4  km  using  the  IASPEI91  travel  times. 
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Figure  3-14.  Comparison  of  regional  seismic  event  location  accuracies  for  14  GT0  PNEs 
from  the  FSU  based  on  observations  from  4-6  IMS  or  surrogate  IMS  stations  using  the 
IASPEI91  travel  times  and  using  travel  times  from  the  regional  3-D  velocity  model. 

We  further  tested  the  location  accuracy  for  our  3-D  regional  velocity  model  in  the 
northern  portion  of  the  DSS  region  of  the  FSU  using  observations  from  a  sample  of 
approximately  20  northern  NZ  explosions,  which  were  observed  at  four  Group  1  regional 
stations,  as  shown  in  Figure  3-15.  The  3-D  DSS  velocity  model  results  again  produce 
big  reductions  in  the  average  mislocations  (15.8  km  versus  43.8  km),  as  illustrated  by  the 
mislocation  vectors  in  Figure  3-16.  As  seen  in  this  figure,  the  event  locations  have  a 
systematic  bias,  because  of  using  only  Group  1  stations,  which  are  all  located  east  of  NZ. 
Adding  data  from  a  single  Group  2  IMS  station  (viz.  OBN),  located  to  the  west,  reduces 
the  bias  and  provides  average  location  accuracy  of  better  than  10  km,  as  can  be  seen  in 
Figure  3-17. 
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Figure  3-15.  Map  locations  of  Group  1  area  regional  stations  (triangles)  which  have 
reported  initial  P  wave  travel  times  from  various  underground  nuclear  tests 
conducted  at  the  Russian  northern  Novaya  Zemlya  (NNZ)  test  site.  A  sample  of 
NNZ  explosions  recorded  at  3  or  more  of  these  stations  has  been  used  for  regional 
seismic  location  analysis  to  further  test  the  applicability  of  our  preliminary  3-D 
velocity  model  in  this  portion  of  the  Group  1  study  area. 
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Figure  3-16.  Comparison  of  mislocation  vectors  for  selected  northern  Novaya  Zemlya 
explosions  computed  with  respect  to  seismic  locations  determined  with  (right)  and  without 
(left)  corrections  to  the  observed  regional  travel  times  with  respect  to  IASPEI91  predicted 
by  the  preliminary  3-D  velocity  model  of  the  Group  1  study  area.  The  arrows  point  from  the 
JED  locations  of  Lilwall  and  Marshall  (1986)  to  the  corresponding  seismic  locations.  Note 
that  the  use  of  the  3-D  velocity  model  results  in.a  very  significant  reduction  in  the  average 
mislocations  (A)  relative  to  IASPEI91  and  that  both  sets  of  mislocation  vectors  are  oriented 
predominantly  west/southwest,  consistent  with  the  uneven  azimuthal  distribution  of  the 
regional  stations  employed  in  the  relocation  analysis. 
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Figure  3-17.  Addition  of  data  from  a  single  regional  station  (OBN)  located  in  the 
Group  2  region  (top)  to  the  west  of  NZ  produces  a  further  significant  reduction  in  the 
average  mislocation  (  A )  of  these  NZ  explosions  (bottom)  to  a  value  of  less  than  10  km. 


So,  the  results  presented  here  for  our  3-D  regional  velocity  model  from  the  FSU  show 
significantly  better  seismic  location  estimates  than  for  the  IASPEI91  nominal  model.  We 
regard  this  as  validation  of  the  model  and  the  associated  regional  travel  time  predictions 
for  the  stations  from  the  DSS  region  of  the  FSU.  Therefore,  we  are  confident  that  the 
SSSCs  predicted  for  the  Group  1  IMS  stations  within  this  region  for  the  3-D  regional 
velocity  model  will  enable  more  accurate  event  locations. 

3.2  CHINA  VALIDATION 

As  described  in  Section  2  above,  our  final  velocity  model  for  the  China/Mongolia  region 
was  developed  from  a  formal  tomographic  inversion  of  observed  travel  times  from  a 
large  sample  of  events  in  the  Annual  Bulletin  of  Chinese  Earthquakes.  We  have 
conducted  a  number  of  analyses  to  validate  this  refined  China/Mongolia  velocity  model. 
One  powerful  means  for  assessing  accuracy  of  the  various  velocity  models  for  China  is 
provided  by  the  recently  published  GT  locations  of  selected  Lop  Nor  explosions  reported 
by  Fisk  (2002).  Figure  3-18  shows  a  comparison  of  Chinese  Bulletin  locations  with  the 
GT  locations  reported  by  Fisk  et  al.  It  can  be  seen  that  the  average  mislocation  is  about 

7.3  km,  with  a  bias  dominantly  to  the  north,  which  reflects  sparse  station  coverage  to  the 
south  of  the  Lop  Nor  test  site.  Although  the  Lop  Nor  epicenters  are  very  well  constrained 


43 


10  km 


* 


88.4*  88.5*  88.6“  88.7*  88.8* 


Figure  3-18.  Mislocation  of  Chinese  Bulletin  epicenters  with  respect  ground 
truth  (GT1-GT2)  locations  for  selected  Lop  Nor  nuclear  explosions. 

by  seismic  and  satellite  image  data,  it  is  necessary  to  assess  the  accuracy  of  the 
corresponding  seismic  origin  times  in  order  to  use  the  seismic  data  from  the  these 
explosions  to  definitively  test  various  velocity  models.  Consequently,  we  conducted  an 
analysis  of  the  reported  CSB  bulletin  arrival  times  of  these  Lop  Nor  explosions  from 
seven  reporting  regional  stations  within  6°.  The  stations  and  their  locations  with  respect 
to  the  Lop  Nor  test  site  are  illustrated  in  Figure  3-19. 

Figure  3-20  shows  a  typical  result  comparing  the  P-wave  travel  time  residuals,  relative  to 
the  origin  times  reported  by  Fisk  et  al.,  for  a  Lop  Nor  explosion  of  06/08/96.  The 
residuals  on  the  left  are  for  the  IASPEI91  travel  times  compared  to  the  residuals  on  the 
right  calculated  from  our  3-D  velocity  model  for  the  region.  In  these  plots  the  symbol 
size  is  proportional  to  the  size  of  the  travel-time  residuals,  and  they  are  plotted  at  the 
epicenter  locations  of  the  corresponding  events.  We  find  a  consistent  bias  of  from  1  to  2 
seconds  at  these  stations  for  the  IASPEI91  model,  but  the  residuals  for  our  3-D  model  are 
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Figure  3-19.  CSB  stations  close  to  Lop  Nor  (A  <  6°)  with  travel  times 
reported  in  the  CSB  bulletin  for  Lop  Nor  explosions. 
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Figure  3-20.  P-wave  travel-time  residuals  at  close  Chinese  stations  for  June  8,  1996  Lop 
Nor  explosion  for  IASPEI91  travel  times  (left)  and  our  3-D  regional  velocity  model  (right). 
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close  to  zero.  More  complete  results  for  12  Lop  Nor  explosions  are  shown  in  Figure  3- 
2 1 .  This  comparison  shows  that  the  average  P-wave  travel  time  residuals  for  these  Lop 
Nor  explosions  with  respect  to  IASPEI91  travel  times  is  about  +  1.25  seconds,  but  the 
corresponding  average  residual  for  our  3-D  model  is  less  than  -  0.2  seconds.  The  scatter 
for  both  sets  of  observations  is  about  the  same  with  a  «  0.3  seconds.  We  conclude  that 
the  Fisk  et  al.  origin  times  for  the  Lop  Nor  events  are  probably  good  to  about  ±  0.5 
seconds.  So,  our  3-D  model  appears  to  account  very  well  for  the  observed  P-wave  travel 
times  in  this  region  close  to  Lop  Nor. 

We  decided  to  look  more  closely  at  the  model  validation  for  the  China  region  using  the 
observations  from  Lop  Nor  explosions  recorded  at  35  regional  stations  with  travel  times 
reported  in  the  CSB  for  multiple  explosions.  The  P-wave  residuals  for  these  events  are 
plotted  in  Figure  3-22.  Although  the  predictions  from  the  3-D  model  look  very  good  out 
to  about  12°,  the  results  appear  to  significantly  deteriorate  at  greater  distances.  At  the 
larger  distances  the  predicted  travel  times  are  generally  too  large,  so  the  model  velocities 
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Figure  3-21 .  Comparison  of  average  P-wave  travel  time  residuals  for  12  Lop  Nor 
explosions  for  IASPEI91  (solid  symbols)  and  for  regional  velocity  model  (open  symbols). 
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Figure  3-22.  Average  Lop  Nor  travel  time  residuals  at  CSB  stations  as  a 

function  of  distance  with  respect  to  IASPEI91  and  3-D  velocity  models. 

affecting  this  distance  range  are  apparently  too  slow.  This  same  behavior  appears  to  be 
corroborated  by  regional  earthquake  travel  time  residuals  from  the  CSB  bulletin  (cf. 
Figure  3-23).  In  this  plot  we  have  binned  the  data  into  0.5°-distance  intervals  and 
computed  the  average  P-wave  residuals.  We  find  a  very  similar  trend  to  the  Lop  Nor 
explosion  data,  again  apparently  indicating  that  the  upper  mantle  velocities  in  the  model 
(below  the  lid)  may  be  too  low. 

We  have  subsequently  revised  the  China/Mongolia  upper  mantle  velocities  in  our  model 
to  help  compensate  for  this  deficiency.  In  particular,  we  analyzed  the  travel  time  data 
from  the  events  in  the  vicinity  of  Lop  Nor  and  determined  that  significantly  higher 
velocities  would  be  required  to  match  the  apparent  Pn  velocities  out  to  larger  regional 
distances.  We  found  that  this  problem  was  apparently  introduced  by  the  relatively  low 
sub-lid  velocities  in  our  base  model  for  the  region  (derived  originally  from  surface-wave 
analyses)  being  carried  through  the  inversion  process.  To  correct  the  problem  we  reran 
the  tomographic  inversion  starting  with  more  uniform,  higher  velocities  in  the  upper 
mantle  for  the  China/Mongolia  base  model.  In  particular,  for  our  base  model  we  imposed 
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Figure  3-23.  Average  CSB  earthquake  travel-time  residuals  as  a  function 
of  distance  with  respect  to  our  China  regional  3-D  velocity  model. 


a  velocity  of  8.12  km/sec  (as  indicated  by  the  data  around  Lop  Nor)  for  the  upper  mantle 
lid  (50-  120  km),  underlain  by  a  uniform  velocity  increase  from  8.12  to  8.3  km/sec 
between  120  and  165  km,  velocity  increasing  from  8.3  to  8.43  km/sec  between  165  and 
210  km,  and  transitioning  to  the  IASPEI91  P  velocity  at  the  410  km  upper  mantle 
discontinuity.  We  reran  the  tomographic  inversion  using  the  same  Chinese  earthquake 
bulletin  data  described  above  with  this  new  starting  base  model  and  determined  a  new 
refined  model  for  China/Mongolia.  The  new  model  has  very  similar  P  velocities  for  the 
upper  mantle  lid  (i.e.  Pn)  as  we  saw  above  in  Figure  2-9.  However,  the  new  model  has 
significantly  higher  sub-lid  P  velocities  throughout  most  of  the  region.  As  a  check  on  this 
model,  we  recomputed  the  P-wave  travel  time  residuals  for  the  Lop  Nor  explosions  at  the 
same  35  stations  shown  in  Figure  3-22.  These  results  are  shown  in  Figure  3-24  based  on 
the  predictions  from  the  new  refined  China/Mongolia  model.  The  predictions  in  this  case 
clearly  match  the  observations  much  better,  and  the  tendency  toward  longer  predicted 
times  at  larger  regional  distances  (A  >  12°)  is  greatly  reduced.  So,  we  have  concluded 
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Figure  3-24.  Average  Lop  Nor  travel  time  residuals  at  CSB  stations  as  a  function 
of  distance  with  respect  to  IASPEI91  and  the  revised  3-D  velocity  model. 

that  the  refined  China/Mongolia  model  developed  from  this  later  tomographic  inversion 
better  represents  the  observations;  and  we  have  substituted  this  new  version  into  our 
overall  3-D  velocity  model  for  the  Group  1  region. 

3.3  WINPAK  VALIDATION 

There  is  very  little  ground  truth  data  currently  available  in  the  India-Pakistan  region  that 
can  be  used  to  validate  our  WINPAK3D  model.  However,  on  14  February  1977,  a 
magnitude  5.2  earthquake  occurred  in  the  region  near  Nilore,  Pakistan  that  was  well 
recorded  by  the  Tarbela  and  Chasma  local  networks.  Based  on  data  from  these  networks, 
Seeber  and  Armbruster  (SA)  (1979)  found  a  hypocenter  for  the  event  and  conducted  a 
detailed  study  of  the  aftershock  sequence.  Using  teleseismic  as  well  as  regional  data, 
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both  the  ISC  and  the  USGS  (NEIC)  located  the  event  within  about  5  km  of  the  epicenter 
reported  by  SA,  who  used  only  the  local  network  data.  The  ISC  calculated  the  depth  of 
the  event  to  be  27  km,  and  the  USGS  fixed  the  depth  at  33  km.  However,  the 
hypocenters  of  the  main  shock  (14.5  km)  and  the  50  accurately  located  aftershocks 
determined  by  S A  indicate  a  rupture  surface  between  1 2  and  1 8  km  depth. 

Because  of  the  regional  station  coverage  and  the  further  constraint  on  the  hypocentral 
region  derived  from  the  aftershock  study,  the  14  February  1977  event  (denoted  the 
VDAY  event)  has  one  of  the  best  constrained  locations  in  the  region.  Therefore,  we 
began  preliminary  testing  of  the  location  capability  of  our  updated  velocity  model  using  a 
subset  of  regional  stations  from  this  event.  Figure  3-25a  shows  the  regional  stations  that 
were  used  to  locate  the  VDAY  event;  the  maximum  azimuthal  gap  is  140°.  We  used  3-D 
travel  time  tables  from  both  our  initial  and  final  WINPAK3D  models  to  locate  the  VDAY 
event  and  compared  the  results  with  the  location  using  the  1-D  IASPEI91  tables.  Figures 
3-25b  and  3-25c  show  the  hypocenters  from  the  VDAY  event,  calculated  using  data  from 
1 1  regional  stations,  for  the  updated  WINPAK3D  model  (blue  *),  initial  WINPAK3D 
model  (red  *),  and  the  1-D  IASPEI91  model  (green  *).  We  compare  these  solutions  with 
the  SA  location  for  this  event  (black  star).  The  surface  and  depth  projections  of  the 
respective  three-dimensional  confidence  regions  determined  by  Monte  Carlo  simulation 
(Rodi  and  Toksoz,  2000)  show  the  95%  confidence  levels  for  each  model's  hypocenter. 
The  epicenter  mislocation  of  the  updated  WINPAK3D  model  from  the  SA  event  location 
is  3.5  km,  while  the  initial  model's  epicenter  mislocation  is  9.3  km  and  the  IASPEI91 
epicenter  mislocation  is  3 1 .7  km.  In  addition,  the  95%  confidence  regions  for  both  the 
initial  and  updated  WINPAK3D  model  epicenters  encompass  the  SA  epicenter  while  the 
95%  confidence  region  for  the  LASPEI91  epicenter  does  not. 

Also  shown  in  Figures  3-25b  and  3-25c  are  the  hypocenters  determined  by  both  the  ISC 
and  the  USGS  (NEIC)  (open  circles)  using  teleseismic  data  as  well  as  regional  data.  Note 
that  by  using  regional  data  alone,  the  3-D  model  is  able  to  estimate  the  hypocenter  of  this 
event  as  well  as  the  teleseismically-derived  estimates.  The  results  of  this  focused  study 
for  an  individual  event  show  that  our  WINPAK3D  regional  velocity  model  can  provide 
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Figure  3-25.  (a)  Eleven  regional  stations  used  to  locate  the  Valentine’s  Day  event. 
Hypocenter  solutions  calculated  with  the  updated  WINPAK3D  model  (blue  *),  the 
initial  WTNPAK3D  model  (red  *),  and  the  1-D  IASPEI91  model  (green  *)  are 
shown  in  b)  map  view  and  c)  depth  vs.  latitude  plots.  The  Seeber  and  Armbruster 
location  for  this  event  is  denoted  by  the  black  star.  Also  shown  are  the  hypocenters 
determined  by  both  the  ISC  and  the  USGS  (NEIC)  (open  circles)  using  teleseismic 
data  as  well  as  regional  data. 

more  accurate  locations  of  small  events,  that  are  not  recorded  teleseismically.  To  further 
investigate  these  capabilities  and  validate  our  model,  we  performed  a  series  of  location 
analyses  on  a  set  of  ground-truth  events.  Our  multi-event  validation  tests  were  performed 
using  a  set  of  17  reference  events  shown  in  Figure  3-26  with  GT5  accuracy  or  better. 
Almost  all  of  the  earthquake  reference  event  locations  were  selected  from  both  the 
Engdahl  and  Bergman  (2001)  cluster  analysis  database  and  from  the  earlier  data  sets 
resulting  from  studies  by  Engdahl  et  al.  (1998).  There  are  three  earthquake  clusters  in  the 
validation  data  set:  Koyna,  Garm,  and  Chamoli.  In  addition  to  the  earthquake  clusters, 
we  include  the  Valentine’s  Day  event,  discussed  above,  and  the  Bhuj,  India  mainshock. 
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Figure  3-26.  Distribution  of  the  set  of  GT5  reference  events  used  for  validation 
of  the  WINPAK3D  model. 

We  also  include  some  explosion  events  in  our  validation  data  set.  These  are  limited  to 
the  two  Indian  nuclear  explosions,  one  Pakistani  nuclear  explosion,  and  two  Peaceful 
Nuclear  Explosions  located  in  the  Former  Soviet  Union,  all  of  which  are  considered  to  be 
known  to  within  1  km.  The  reference  event  data  set  was  intentionally  excluded  from  our 
inversion  database  to  allow  for  an  independent  validation  analysis. 

Our  model  validation  approach  included  relocation  of  our  set  of  reference  events  and 
analysis  of  the  results  for  three  cases:  fixed  hypocenter,  fixed  depth,  and  no  constraints. 
We  present  the  results  in  Figures  3-27  through  3-29,  identifying  the  number  of  events 
improved  versus  the  number  of  events  deteriorated  relative  to  IASPEI91 .  The  graphs 
show  percentage  improvement  (epicenter,  depth,  RMS  travel  time  residuals)  over 
IASPEI91  versus  number  of  events.  When  we  fixed  the  hypocenters  to  the  GT5 
locations  (see  Figure  3-27),  we  found  that  while  the  mean  travel  time  residuals  for  both 
the  IASPEI91  model  and  the  WINPAK3D  model  were  approximately  zero,  the  standard 
deviation  was  1.85  sec  for  the  WINPAK3D  model  and  2.48  sec  for  the  IASPEI91  model. 
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Figure  3-27.  Validation  results  for  the  case  where  reference  events  are  fixed  to 
their  known  locations:  (a)  Travel  time  residuals  for  IASPEI91  (green  diamonds) 
and  the  WINPAK3D  (blue  dots)  models,  (b)  Percentage  improvement  in  RMS 
travel  time  error  compared  with  IASPEI91  versus  number  of  events. 
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Figure  3-28.  Validation  results  for  the  case  where  epicentral  location  is  determined, 
but  the  depths  of  reference  events  are  fixed.  Percentage  improvement  compared  to 
IASPEI91  for  (a)  epicentral  mislocation  and  (b)  RMS  travel  time  error. 
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Figure  3-29.  Validation  results  for  free  hypocenter  relocations.  Percentage 
improvement  compared  to  IASPEI91  in  (a)  epicentral  mislocation,  (b)  depth 
mislocation,  and  (c)  RMS  travel  time  error. 


In  general,  the  IASPEI91  residuals  exhibit  a  greater  bias  towards  positive  residuals  out  to 
10  degrees  in  this  region,  with  a  negative  bias  beyond  10  degrees.  Additionally,  the 
WINPAK3D  model  reduced  RMS  travel  time  error  for  all  17  events  by  a  mean  of  21.2%, 

In  the  second  case,  we  fixed  the  event  depths  and  solved  for  epicentral  locations.  The 
results  are  shown  in  Figure  3-28.  While  WINPAK3D  improved  the  RMS  travel  time 
residuals  for  16  out  of  17  events,  it  improved  epicentral  mislocation  for  only  12  out  of  17 
events,  or  70%  of  the  events.  Mean  percent  improvement  in  epicentral  mislocation  was 
33.6%.  For  those  events  where  event  locations  were  degraded,  IASPEI91  mean  percent 
improvement  was  38.7%.  With  regard  to  RMS  travel  time  residuals,  mean  percent 
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improvement  was  23.7%  for  WINPAK3D.  For  the  one  event  that  had  a  lower  RMS  error 
based  on  the  IASPEI91  model,  its  mean  percent  improvement  over  the  WINPAK3D 
model  was  only  9.7%. 

In  the  last  case,  we  performed  a  free  hypocenter  relocation.  These  results  are  shown  in 
Figure  3-29.  Epicentral  mislocation  was  improved  for  12  of  17  events  with  a  mean 
improvement  of  38.7%.  For  the  5  events  where  the  IASPEI91  model  found  more 
accurate  locations,  its  mean  improvement  over  the  WINPAK3D  model  was  39.4%. 

Depth  mislocation,  in  this  case,  demonstrates  the  most  dramatic  improvement.  The 
WINPAK3D  model  improved  depth  accuracy  for  13  of  17  events  by  a  mean 
improvement  of  68.5%.  On  the  other  hand,  the  mean  improvement  was  only  18.0%  for 
the  4  events  where  IASPEI91  found  a  more  accurate  event  depth.  Finally,  RMS  travel 
time  error  was  improved  for  14  of  17  events  by  a  mean  of  24.5%.  The  IASPEI91  model 
had  smaller  RMS  residuals  for  3  events  with  a  mean  improvement  over  WINPAK3D  of 
30.8%. 


These  results  demonstrate  that  the  WINPAK3D  model  is  a  significant  improvement  over 
the  LASPEI91  model  for  regional  location  calibration.  For  this  set  of  17  reference  events, 
travel  times  calculated  through  the  WINPAK3D  model  were  more  accurate  than  those 
calculated  with  the  IASPEI91  model  for  at  least  80%  of  the  data.  As  a  result,  event 
locations  determined  using  the  WINPAK3D  model  were  found  to  be  more  accurate,  in 
both  depth  and  epicenter,  for  at  least  70%  of  the  data  set.  We  are  confident  that  the 
WINPAK3D  sub-model  is  a  substantial  improvement  in  regional  calibration  for  the 
Group  1  IMS  stations  in  and  around  the  India-Pakistan  region. 
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SECTION  4 


EXTENSION  TO  SECONDARY  SEISMIC  PHASES 
4.1  CONVERSION  TO  S-WAVE  VELOCITY  MODELS 

In  addition  to  the  primary  (P)  phases,  secondary  seismic  phases  may  be  especially  useful 
in  locating  small  regional  seismic  events.  For  example,  traditional  methods  of  locating 
regional  events  often  utilize  S  -  P  times  to  provide  preliminary  estimates  of  source-to- 
station  distance.  More  sophisticated  seismic  location  algorithms  can  use  a  variety  of 
secondary  seismic  phase  times  to  achieve  more  reliable  location  estimates.  However,  just 
as  with  primary  phases,  the  secondary  phase  travel  times  need  to  be  calibrated  on  a 
regional  basis  to  be  helpful  in  reducing  location  error.  To  help  in  regional  calibration  of 
secondary  seismic  phases,  we  have  conducted  preliminary  analyses  of  model-based 
methods,  following  procedures  similar  to  those  described  above  for  the  regional  P  waves. 
In  particular,  to  test  this  approach  we  have  derived  and  conducted  initial  validation  of  a  3- 
D  S-wave  velocity  model  for  the  Group  1  region  by  converting  from  the  P-wave  model 
using  relationships  based  on  Poisson’s  ratio  values.  This  3-D  S-velocity  model  can  then 
be  used  to  predict  regional  travel  times  and  corrections  (SSSCs)  for  S  and  possibly  other 
secondary  phases  (e.g.  Lg),  which  may  be  useful  for  event  location. 

We  assume  that  the  S  velocity  model  has  the  same  general  description  as  our  3-D  P 
velocity  model  represented  by  the  UMP  velocity-versus-depth  profiles  at  the  geographic 
grid  elements  spaced  uniformly  in  latitude  and  longitude.  The  S  velocities  at  each  node 
in  the  model  are  simply  obtained  from  the  relationship: 


VS=VP 


l-2o 

2(1-°) 


1/ 
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where  a  is  the  Poisson’s  ratio  and  Vp  the  P  velocity  at  any  node  in  the  model. 


(4-1) 


4.2  POISSON’S  RATIOS  FOR  THE  CRUST  AND  UPPER  MANTLE 
Values  of  Poisson’s  ratio  for  rocks  in  the  earth’s  crust  and  upper  mantle  are  generally  in 
the  range  of  a  «  0.25.  Rocks  in  the  earth’s  crust  generally  have  somewhat  lower 
Poisson’s  ratios  than  rocks  in  the  upper  mantle.  Thus,  for  the  IASPEI91  model  Poisson’s 
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ratios  fall  in  the  range  from  about  0.24  to  0.29  over  depths  from  the  earth’s  surface  down 
to  the  upper  mantle  velocity  discontinuity  near  410  km.  Although  this  range  does  not 
appear  to  be  very  large,  it  should  be  noted  that  small  differences  in  the  Poisson’s  ratios 
can  produce  large  differences  in  the  corresponding  S  velocities  within  the  earth  model 
and  in  the  associated  S-wave  travel  time  calculations. 

Toksoz  and  his  colleagues  at  MIT  reviewed  published  reports  for  the  existing  P-  and  S- 
wave  velocity  models  for  China  and  analyzed  the  relationship  between  P  and  S  velocities 
for  the  crust  and  upper  mantle.  In  general,  they  found  differences  in  the  velocity  ratios 
between  the  crust  and  upper  mantle.  In  particular,  they  found  that  Vp/Vs  ratios  in  China 
increase  with  depth  and  were  in  the  range  from  1.74  to  1.83  at  crustal  depths  and  in  the 
range  from  1.83  to  1.85  in  the  upper  mantle.  The  corresponding  Poisson’s  ratios  have 
values  near  0.21  for  crustal  rocks  with  P  velocities  near  5.50  km/sec  and  increase  to 
values  near  0.29  for  crustal  and  upper  mantle  rocks  with  P  velocities  near  7.80  km/sec 
and  higher.  The  MIT  group  results  indicate  that  the  relationship  between  the  P  and  S 
velocities  for  China  crust  and  upper  mantle  rocks  can  be  expressed  in  the  form  of  a 
simple  relationship: 

=  2.04  VP -1.655  (4-2) 

So,  within  China  we  use  this  relationship  to  derive  the  S  velocities  from  the  P  velocities 
in  our  3-D  model. 

For  the  model  outside  of  China,  we  adopted  a  constant  Poisson’s  ratio  of  a  =  0.26,  which 
implies  Vp/Vs  =  1 .7559  everywhere.  Under  this  assumption  the  S-wave  raypaths  through 
the  3-D  model  are  exactly  the  same  as  the  corresponding  P-wave  raypaths,  and  Ts/TP  = 

1 .7559.  We  can  use  this  to  derive  a  simple  relationship  between  the  S-wave  SSSCs  and 
the  P-wave  SSSCs,  as  follows.  By  definition 


SSSC(P)  = 

Tp  (3-DModel)-TP  (IASPEI91) 

(4-3) 

sssc(s)  = 

Ts  (3-DModel)-Ts  (IASPEI91) 

(4-4) 

= 

1.7559  TP  (3-DModel)-Ts(lASPEI9l) 

(4-5) 
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1.7559 


Ts  (IASPEI91) 


(4-6) 


TP  (3-D  Model) 


1.7559 


=  1.7559 


SSSC(P)+TP  (IASPEI91)- 


Ts  (IASPEI91 ) 
1.7559 


(4-7) 


So,  for  this  constant  Poisson’s  ratio,  the  S-wave  SSSCs  can  be  expressed  simply  in  terms 
of  the  P-wave  SSSCs  as: 


SSSC(S)  =  1.7559  SSSC(p)+  [l.7559TP  (lASPEI9l)-Ts  (IASPEI91 )]  (4-8) 

For  areas  where  this  constant  Poisson’s  ratio  is  applicable,  we  can  then  use  the  calculated 
P-wave  SSSCs  to  generate  the  corresponding  S-wave  SSSCs  by  applying  the  simple 
correction  involving  the  IASPEI91  P  and  S  travel  times. 

We  illustrate  this  S-wave  travel  time  prediction  capability  for  a  path  in  the  FSU  platform 
region  to  the  northeast  of  station  Borovoye  in  Figure  4-1 .  This  plot  compares  the  S-wave 


SSSC(S)  from  3-D  Model 

Figure  4-1 .  Comparison  of  the  prediction  of  S-wave  SSSCs  for  the  region  northeast 
of  Borovoye  based  on  the  relationship  of  Equation  4-8  for  a  constant  Poisson’s  ratio 
with  those  determined  by  raytracing  through  the  3-D  velocity  model. 


58 


SSSCs  predicted  by  the  relationship  in  Equation  4-8  with  the  S-wave  SSSCs  calculated 
using  the  raytracing  algorithm  through  the  3-D  S-wave  velocity  model  along  the  45° 
azimuth  from  station  BRVK.  The  plot  shows  the  SSSC  predictions  from  Equation  4-8  on 
the  vertical  axis  compared  to  the  model-based  SSSCs  on  the  horizontal  axis  over  the 
distance  range  from  3°  <  A  <  10°.  The  comparison  shows  that  the  predicted  SSSCs  using 
the  simple  relationship  based  on  a  constant  Poisson’s  ratio  match  almost  exactly  the  times 
predicted  by  raytracing  through  the  3-D  model  over  this  distance  range  for  this  region. 

So,  in  regions  where  we  can  make  this  simple  assumption  on  Poisson’s  ratio,  we  have  a 
very  simple  alternative  to  predict  the  S-wave  SSSCs.  In  particular,  we  do  not  need  to  do 
the  additional  raytracing  for  the  S-waves;  but  instead  we  can  use  the  P-wave  SSSCs  with 
the  modification  from  Equation  4-8. 

4.3  PRELIMINARY  VALIDATION  OF  S-WAVE  MODELS 
We  have  conducted  a  number  of  preliminary  analyses  and  tests  to  demonstrate  the 
reasonableness  and  validity  of  the  3-D  S-wave  velocity  models  for  several  areas  within 
the  Group  1  study  region.  One  such  validation  test  was  for  the  region  around  station 
Borovoye  (BRVK).  Figure  4-2a  illustrates  two  great  circle  paths  extending  out  to  20°  at 
azimuths  due  east  and  northeast  (45°)  from  BRVK.  We  used  the  PL  raytracing  algorithm 
and  our  3-D  S-wave  velocity  model  to  calculate  the  S-wave  travel  times  along  these 
profiles.  The  results  are  shown  in  Figure  4-2b.  We  see  some  differences  in  the  travel 
times  along  the  two  profiles  with  apparent  average  S  velocities  (estimated  by  linear  fits  to 
the  calculated  travel  times)  of  4.44  km/sec  along  the  due  east  profile  and  4.66  km/sec 
along  the  northeast  profile.  In  general,  such  differences  look  reasonable  considering  that 
we  would  expect  somewhat  slower  S  velocities  along  the  east  profile  which  runs  partially 
through  tectonically  active  areas. 

In  addition,  we  plotted  S-wave  travel  times  assembled  by  the  IDG  from  FSU  bulletin  data 
for  PNEs  recorded  at  various  stations  of  the  Soviet  permanent  network.  These  are  plotted 
in  Figure  4-3  over  the  regional  distance  range  out  to  20°.  We  again  fit  a  simple  linear 
straight  line  relationship  to  these  data  and  obtained  an  apparent  Sn  velocity  of  4.65 
km/sec.  This  is  almost  exactly  equal  to  the  apparent  S  velocity  seen  above  in  our  model 
for  the  profile  northeast  from  BRVK.  So,  the  regional  S-wave  observations  from  the 
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Figure  4-2a.  Great  circle  paths  extending  out  to  20°  to  the  east  (red)  and  to  the 
northeast  (blue)  from  the  Borovoye  station  used  for  computing  the  S-wave  travel 
times  from  the  3-D  regional  model. 


A, degrees 


j.  iS^avertravehtimexurvesxomputexhusmg^e-PB-algorilhm-for 
profiles  to  the  east  (red)  and  northeast  (blue)  from  the  Borovoye  station  through 
the  3-D  velocity  model. 
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Figure  4-3.  Observed  Sn  travel  times  versus  distance  for  events  recorded  in 
the  FSU  from  GTO  PNEs. 

FSU  PNEs  appear  to  generally  corroborate  the  Sn  velocities  predicted  by  our  model  from 
the  region  around  BRVK. 

In  Figure  4-4  we  show  a  more  detailed  comparison  for  the  S-waves  observed  at  Borovoye 
from  FSU  PNEs.  We  compare  here  the  S-wave  residuals  with  respect  to  IASPEI91  (on 
the  left)  and  with  respect  to  our  3-D  S  velocity  model  (on  the  right).  We  see  that  there  is 

a  very  large  negative  bias  ( At  =  -  8.83  sec )  for  the  observed  S  residuals  relative  to  the 

IASPEI  travel  times.  The  bias  is  essentially  removed  ( At  =  1 .54  sec)  by  the  3-D  S 
velocity  model  results  seen  in  the  right-hand  figure.  Although  these  data  are  limited, 
there  do  appear  to  be  some  systematic  trends  in  the  residuals  with  respect  to  the  3-D 
model  which  suggest  that  a  more  complete  analysis,  including  a  tomographic  inversion  to 
recover  a  spatially  dependent  Poisson’s  ratio,  might  be  useful. 
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Figure  4-4.  Comparison  of  S-wave  travel  time  residuals  at  station  Borovoye  for 
PNEs  measured  with  respect  to  IASPEI91  travel  times  (left)  and  with  respect  to  our 
3-D  S-wave  velocity  model. 

One  thing  that  does  seem  clear  at  this  point  is  that  the  variance  in  the  S-wave  arrival 
times  will  be  significantly  larger  than  that  for  P,  even  for  a  very  good  S-wave  velocity 
model.  For  example,  for  this  Borovoye  sample,  if  we  estimate  the  best-fitting  average 
value  of  Poisson’s  ratio  by  averaging  the  observed  ratios  of  S-to-P  travel  times  and  use 
this  ratio  to  predict  S-wave  travel  times  from  the  observed  P-wave  travel  times  (as 
described  above),  the  standard  deviation  of  the  difference  between  the  predicted  and 
observed  S  times  is  about  4  seconds,  as  opposed  to  about  1  second  for  the  P-wave  times 
predicted  by  our  3-D  P-wave  velocity  model.  This  is  the  best  that  can  be  expected  from  a 
constant  Poisson’s  ratio  conversion  from  our  P-wave  model,  and  it  suggests  that  S-wave 
arrival  times  should  be  significantly  downweighted  in  the  location  process,  at  least  for 
explosions. 

A  related  problem  common  to  all  secondary  phases  is  the  ambiguity  in  determining  onset 
times.  This  fact  is  illustrated  in  Figure  4-5,  which  shows  a  comparison  of  S-wave  onset 
PNE_evenlsj£c^ 
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Figure  4-5.  Comparison  of  differences  in  S-wave  onset  time  estimates  at 
station  Borovoye  by  analysts  from  Lamont  and  from  IDG  for  FSU  PNEs. 

Lamont  and  at  the  Russian  IDG.  Despite  the  fact  that  these  are  very  senior  analysts 
working  with  high-quality  digital  data,  it  can  be  seen  that  there  are  very  significant 
differences,  averaging  more  than  6  seconds,  between  the  two  sets  of  picks.  This  suggests 
that  secondary  phase  travel  time  corrections  will  have  to  be  very  carefully  calibrated  to 
automatic  processing  and  analyst  procedures  at  the  processing  center  of  interest  (e.g. 
IDG). 

4.4  OTHER  SECONDARY  PHASES 

As  noted  above,  other  secondary  regional  phases,  including  Pg  and  Lg,  can  also  be  useful 
for  determining  better  seismic  locations  for  small  events,  provided  they  are  properly 
calibrated.  For  example,  Figure  4-6  shows  Pg  and  Lg  travel  times  from  FSU  PNEs  which 
would  be  useful  in  calibrating  regional  travel  times  for  these  phases  at  IMS  stations  from 
events  in  our  Group  1  study  region.  As  part  of  this  project,  we  investigated  the 
possibility  of  defining  the  travel-time  corrections  for  these  secondary  regional  phases  by 
mapping  the  travel  times  of  the  first  crustal  arrivals  through  our  P  and  S  wave  models 
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Figure  4-6.  Observed  secondary  phase,  Pg  and  Lg,  travel  times 
versus  distance  for  GTO  PNEs  from  the  FSU. 

using  the  PL  algorithm  and  “masking”  out  the  upper  mantle  arrivals.  Unfortunately,  our 
initial  studies  found  that  the  arrivals  isolated  by  this  approach  appear  to  correspond  to 
deep  crustal  refractions,  which  had  group  velocities  inconsistent  with  the  nominal 
observed  Pg  and  Lg  velocities  and,  consequently,  were  not  obviously  applicable  to  the 
definition  of  the  travel  time  corrections  for  these  phases. 

We  also  looked  at  an  alternate  approach  based  on  published  results  from  various  regional 
synthetic  seismogram  studies,  which  have  indicated  that  the  onset  of  the  Pg  and  Lg 
arrival  groups  correspond  to  the  Moho  reflections  PmP  and  SmS,  respectively.  Since  the 
vertical  P  and  S  wave  travel  times  through  the  IASPEI91  crustal  model  imply  average  P- 
and  S-wave  crustal  velocities  of  about  6,1  and  3.5  km/sec,  respectively,  the  predicted 
group  velocities  for  these  phases  are  consistent  with  the  IASPEI91  travel  time  curves  for 
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Pg  and  Lg.  We  believe  that  it  should  be  feasible  to  adapt  the  PL  algorithm  to  the 
estimation  of  PmP  and  SmS  travel  times  through  our  3-D  model  for  use  in  the  definition 
of  SSSC  surfaces  for  secondary  phases,  Pg  and  Lg.  Such  calibration  corrections  could 
make  these  regional  phases  much  more  valuable  for  seismic  location  of  small  events 
detected  by  limited  numbers  of  regional  stations.  However,  additional  research  will  be 
required  before  such  corrections  can  be  reliably  applied  to  observed  Pg  and  Lg  arrival 
time  data. 

It  is  also  noteworthy  that,  as  with  the  Sn  discussion  above,  careful  calibration  will  be 
needed  to  tailor  any  such  Pg  and  Lg  travel  time  corrections  to  specific  applications.  That 
is,  the  corrections  proposed  above  strictly  refer  to  the  theoretical  Pg  and  Lg  onset  times, 
which  may  be  quite  different  from  the  Pg  and  Lg  “arrival  times”  normally  determined  by 
analysts  who  typically  focus  on  later,  higher  amplitude  portions  of  the  envelope  of 
arrivals  which  define  these  phase  observationally.  This  secondary  phase  pick  uncertainty 
is  illustrated  in  Figure  4-6  above  which  shows  considerably  greater  scatter  about  the 
average  lines  than  was  seen  for  the  Sn  (cf.  Figure  4-3),  and  even  more  clearly  compared 
to  the  P-wave  times. 


SECTION  5 


SSSCs  FOR  P  AND  S  AT  GROUP  1  IMS  STATIONS 

5.1  DETERMINATION  OF  SSSCs 

The  defined  deliverables  for  this  project  are  travel  time  corrections  for  each  of  the 
identified  Group  1  IMS  stations,  implemented  in  the  form  of  SSSCs,  where  the  SSSC 
represents  the  difference  (as  a  function  of  source  location)  between  the  travel  times 
predicted  by  the  derived  earth  model  for  the  region  and  those  predicted  by  the  default 
IASPEI91  model.  As  part  of  our  program  efforts  over  the  past  year,  we  have  recomputed 
the  SSSCs  based  on  our  final  3-D  velocity  model  for  all  30  IMS  seismic  stations  from  the 
Group  1  study  region.  We  have  now  computed  and  delivered  digital  versions  of  the 
SSSCs  corresponding  to  both  the  P-  and  the  S-wave  travel  time  corrections  (relative  to 
IASPEI91)  extending  out  to  20°  from  each  station.  These  correction  tables  have  been 
specified  for  surface  focus  plus  10  other  focal  depths  (viz.  5,  10,  20,  30,  40,  50,  70,  100, 
150,  and  200  km). 

Plots  of  the  surface-focus  SSSCs  for  all  30  Group  1  IMS  stations  are  presented  in  the 
appendix  to  this  report.  In  the  subsections  which  follow,  we  describe  some  examples  of 
the  P-  and  S-wave  SSSCs  for  selected  Group  1  IMS  stations.  We  also  look  at  the  depth 
dependence  of  the  SSSCs  and  discuss  the  application  of  these  corrections  for  use  in  event 
location  in  the  Group  1  region. 

5.2  EXAMPLES  OF  P-WAVE  SSSCs  FOR  SELECTED  GROUP  1  IMS  STATIONS 
At  the  conclusion  of  the  second  year  of  this  research  program,  we  prepared  estimates  of 
P-wave  SSSCs  for  all  30  of  the  Group  1  IMS  stations  based  on  the  existing  3-D  velocity 
model  for  the  region.  As  noted  in  the  discussion  above,  there  have  been  a  number  of 
revisions  and  refinements  to  the  velocity  model  throughout  this  region  as  a  result  of 
continuing  calibration  efforts  with  additional  ground  truth  data.  In  particular,  the  new 
model  incorporates  major  changes  to  the  China  region,  based  on  our  formal  tomographic 
analysis  of  the  earthquake  data  from  that  region,  as  well  as  a  number  of  minor  changes  in 
some  other  areas  (e.g.  WINPAK3D  refinements  and  increased  upper  mantle  velocities  in 
some  of  the  continental  areas  to  the  far  east  in  the  FSU,  which  had  been  formerly 
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represented  by  surface-wave  models).  The  final  3-D  regional  velocity  model  for  the 
entire  Group  1  region  of  eastern  Asia  was  assembled  using  the  QUILT  algorithm  to 
compose  the  revised  submodels,  as  described  above.  We  then  used  the  PL  algorithm  to 
compute  the  P-wave  travel  times  through  this  composite  3-D  model  for  distances  out  to 
20°  surrounding  each  of  the  IMS  station  sites  for  each  of  the  assumed  focal  depths.  The 
final  SSSCs  are  then  generated  as  the  difference  between  these  computed  values  and  the 
corresponding  P-wave  travel  times  from  the  IASPEI91  tables. 

Figures  5-1  through  5-5  show  representative  examples  of  the  P-wave  SSSCs  for  surface 
focus  events  at  five  selected  IMS  station  sites  in  the  Group  1  region  of  eastern  Asia.  We 
have  plotted  all  of  the  SSSC  values  on  a  common  color  scale  to  facilitate  comparisons. 
The  first  example  is  for  station  Borovoye  (cf.  Figure  5-1).  The  P-wave  SSSC  values  at 
Borovoye  vary  from  -  9.2  seconds  to  +  2.9  seconds.  As  one  might  expect,  the 
corrections  tend  to  be  small  close  to  the  station;  and  the  absolute  values  show  a  tendency 
to  increase  (i.e.  corrections  become  more  negative)  with  distance.  In  general,  the 
corrections  over  much  of  the  region  surrounding  Borovoye  tend  to  be  negative,  which 
indicates  that  the  3-D  model  travel  times  are  shorter  than  the  corresponding  LASPEI91 
times.  This  means  that  the  velocities  in  our  3-D  model  are  generally  higher  than 
IASPEI91  1-D  model  velocities  for  much  of  this  region.  Some  of  the  largest  negative 
corrections  occur  at  north  and  northwesterly  azimuths  (magenta  colors  in  the  plot).  In 
general,  the  correction  map  for  Borovoye  appears  quite  reasonable  and  represents  the 
expected  behavior,  in  that  the  largest  negative  corrections  are  for  paths  which  cross 
continental  platform  areas  which  tend  to  have  higher  than  normal  seismic  velocities.  In 
contrast,  it  also  appears  reasonable  in  Figure  5-1  that  the  corrections  coming  from  the 
active  tectonic  area  along  the  southern  border  of  the  FSU  with  China/Mongolia  tend  to  be 
near  zero  or  positive,  because  these  path  segments  would  tend  to  have  velocities  lower 
than  normal.  Although  the  SSSC  corrections  are  strictly  only  valid  out  to  20°,  the  plots 
in  this  and  subsequent  figures  extend  to  somewhat  greater  distances  in  the  comers;  and 
we  see  in  these  areas  that  the  travel  time  corrections  at  larger  distances  tend  to  return  to 
more  nominal  values  (i.e.  corrections  closer  to  zero).  This  behavior  also  appears 
consistent  with  what  we  expect  from  the  model,  since  at  larger  distances  the  travel  times 
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Figure  5-2.  Estimated  surface-focus  P-wave  SSSCs  for  the  LZDM  station 
site  range  from  -3.2  seconds  to  +3.3  seconds. 
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Figure  5-1.  Estimated  surface-focus  P-wave  SSSCs  for  the  BRVK  station 
site  range  from  -9.2  seconds  to  +3.0  seconds. 


Time  (s) 


68 


Figure  5-3.  Estimated  surface-focus  P-wave  SSSCs  for  the  NRIK  station  site 
range  from  -8.9  seconds  to  +3.3  seconds. 
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figure  5-4.  Estimated  surface-focus  P-wave  SSSCs  for  the  CMTO  station  site  range 
from  -2.9  seconds  to  +5.8  seconds. 
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Figure  5-5.  Estimated  surface-focus  P-wave  SSSCs  for  the  USK  station  site  range 
from  -4.6  seconds  to  +7.3  seconds. 

are  most  affected  by  deeper  parts  of  the  model,  which  tend  to  return  to  IASPEI91  1-D 
values  for  rays  bottoming  near  the  410-km  discontinuity  (i.e.  below  410  km  our  3-D 
model  matches  the  IASPEI91  model). 

Figure  5-2  shows  a  rather  different  behavior  for  the  P-wave  SSSC  values  from  the  region 
around  the  LZDM  station  site  in  central  China.  The  predicted  P-wave  SSSC  values  for 
LZDM  range  only  from  -  3.2  seconds  to  +  3.3  seconds.  We  again  see  that  values  near 
the  station  are  about  zero  and  the  absolute  corrections  again  generally  increase  with 
distance  over  the  range  out  to  20°.  Largest  negative  corrections  occur  near  the  maximum 
regional  distance.  In  general,  the  corrections  around  LZDM  appear  to  be  rather  uniform 
and  do  not  show  great  differences  on  average  with  respect  to  IASPEI91  (i.e.  correction 
values  tend  to  average  near  zero  over  much  of  the  region).  However,  there  again  appears 
to  be  some  evidence  of  a  swath  across  the  southwestern  border  of  China  where  the 
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correction  values  near  20°  appear  to  be  less  negative.  This  again  would  seem  to  be 
consistent  with  generally  lower  velocities  for  paths  across  this  tectonically  active  zone. 

Figure  5-3  shows  the  surface-focus  P-wave  SSSCs  from  the  region  around  station  Norilsk 
(NRIK)  in  northern  Russia.  In  this  case  the  SSSC  values  range  from  -  8.9  seconds  to  3.3 
seconds.  Values  near  the  station  are  again  near  zero  with  some  positive  values  in  the 
nearer  regional  distance  ranges.  The  largest  negative  corrections  occur  near  the  farther 
limits  of  the  regional  distance  range  to  the  southeast  and  south.  In  general,  most  of  the 
raypaths  from  sources  in  these  areas  tend  to  be  across  the  platform  areas  of  the  FSU;  so 
the  model  velocities  in  these  areas  tend  to  be  higher  than  IASPEI91  in  these  areas 
producing  the  expected  negative  SSSCs.  In  this  case,  the  tendency  for  the  tectonic  areas 
along  the  FSU  southern  border  with  Mongolia  to  have  less  negative  correction  values  is 
not  apparent,  perhaps  because  the  corresponding  zone  of  lower  velocities,  which  we 
expect  in  this  area,  is  too  close  to  the  far  extreme  of  the  paths  from  this  area  and 
producing  only  a  minor  effect  on  the  travel  times.  One  result,  which  gives  some 
confidence  in  these  SSSC  predictions,  can  be  seen  by  comparing  Figure  5-3  with  Figure 
5-1  for  the  path  between  the  two  station  locations.  First,  looking  at  the  SSSC  value  in 
Figure  5-3  for  station  NRIK  from  a  potential  event  at  the  earth’s  surface  near  the  location 
of  station  BRVK,  the  correction  value  is  approximately  -  6  seconds  (i.e.  violet  in  the 
color  scale).  Next,  looking  at  Figure  5-1  for  station  BRVK  from  a  potential  event  near 
the  location  of  station  NRIK,  the  correction  value  is  again  approximately  -  6  seconds. 

So,  the  SSSC  corrections  calculated  from  the  models  for  these  two  reciprocal  paths  are 
about  the  same,  as  we  would  expect  from  the  prediction  scheme. 

Figure  5-4  and  5-5  show  two  additional  examples  of  the  surface-focus  P-wave  SSSCs 
predicted  by  the  model  for  stations  CMTO  in  Thailand  and  USK  in  eastern  Russia.  For 
CMTO  the  SSSCs  are  in  the  range  from  -  2.9  seconds  to  +  5.8  seconds,  and  for  USK  the 
SSSCs  are  in  the  range  from  -  4.6  seconds  to  +  7.3  seconds.  Values  close  to  the  stations 
in  both  cases  are  again  near  zero.  The  corrections  with  the  most  negative  values  for  both 
of  these  stations  are  for  continental  paths  at  the  largest  regional  distances.  In  particular, 
most  negative  values  in  Figure  5-4  occur  in  a  relatively  limited  area  to  the  northeast  of 
station  CMTO  at  17°  to  20°;  and  in  Figure  5-5  the  most  negative  values  occur  in  a  wide 
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arc  to  the  west  of  USK  again  at  regional  distances  from  about  17°  to  20°.  A  significant 
difference  in  the  SSSC  plots  in  Figures  5-4  and  5-5  is  the  occurrence  of  relatively  large 
positive  values  from  events  with  paths  through  the  oceanic  areas  surrounding  these 
stations.  Thus,  for  station  CMTO  we  see  large  positive  SSSC  values  from  locations  in 
oceanic  areas  to  the  southwest  and  to  the  far  east;  and  for  station  USK  there  are  large 
positive  SSSC  values  from  location  in  the  oceanic  areas  to  the  far  southeast.  These  SSSC 
prediction  results  again  would  appear  to  be  consistent  with  what  we  expect  from  our 
general  3-D  model,  which  we  know  to  have  low  velocities  relative  to  IASPEI91  in 
oceanic  areas. 

5.3  EXAMPLES  OF  S-WAVE  SSSCs  FOR  SELECTED  GROUP  1  IMS  STATIONS 
In  Section  4  we  described  the  procedures  used  to  generate  the  3-D  S-wave  velocity  model 
from  the  P-wave  model.  The  final  S-wave  3-D  velocity  model  was  then  used  as  input  to 
the  PL  algorithm  again  to  calculate  the  S-wave  travel  times  out  to  20°  from  each  of  the  30 
IMS  station  sites  in  the  Group  1  area  for  each  of  the  assumed  focal  depths.  The  S-wave 
SSSCs  for  each  station  were  then  computed  as  the  difference  between  the  travel  times 
from  the  3-D  regional  model  and  the  corresponding  S-wave  travel  times  from  the 
IASPEI91  tables. 

In  Figures  5-6  through  5-8,  we  show  some  representative  examples  of  the  S-wave  SSSCs 
for  surface  focus  events  at  three  selected  IMS  station  sites  from  the  Group  1  region  of 
eastern  Asia.  Again  all  of  the  plots  are  on  a  common  color  scale  to  facilitate 
comparisons,  although  the  scale  for  the  S  waves  is  different  from  that  shown  above  for  P 
waves  to  show  the  greater  range  in  the  S-wave  corrections.  The  first  plot  is  for  station 
Borovoye  (cf.  Figure  5-6).  The  S-wave  SSSCs  predicted  by  our  model  have  values  at 
Borovoye  from  -  29.5  seconds  to  +  5.4  seconds.  As  seen  in  the  plot,  the  values  are  close 
to  zero  near  the  station  and  the  absolute  values  generally  increase  with  distance  away 
from  the  station.  The  most  negative  corrections  occur  to  the  north  and  northwest  at  far 
regional  distances.  The  plot  for  the  S-wave  SSSCs  in  Figure  5-6  again  shows  evidence  of 
the  low-velocity  tectonic  zone  along  the  southern  border  of  the  FSU  with 
China/Mongolia  in  the  form  of  a  swath  of  somewhat  more  positive  correction  values. 
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Figure  5-6.  Estimated  surface-focus  S-wave  SSSCs  for  the  BRVK  station  site 
range  from  -29.5  seconds  to  +5.4  seconds. 
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Figure  5-7.  Estimated  surface-focus  S-wave  SSSCs  for  the  LZDM  station  site 
range  from  -25.6  seconds  to  +46.9  seconds. 
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Figure  5-8.  Estimated  surface-focus  S-wave  SSSCs  for  the  CMTO  station 
site  range  from  -24.1  seconds  to  +42.1  seconds. 

The  S-wave  SSSC  plot  in  Figure  5-7  for  station  LZDM  in  central  China  exhibits  much 
the  same  kind  of  uniformity  that  we  saw  above  for  the  P-wave  SSSCs  at  this  station. 
Although  the  predicted  surface-focus  S-wave  SSSCs  for  LZDM  range  from  -  25.6 
seconds  to  +  46.9  seconds,  the  extreme  positive  and  negative  values  all  come  from  the 
comers  of  the  plot  beyond  the  20°  regional  distance  range.  In  fact,  within  the  20°  range 
we  see  that  the  predicted  S-wave  SSSCs  are  quite  moderate  (in  the  range  from  about  -  9 
seconds  to  +  9  seconds)  and  there  are  only  minor  fluctuations. 

Figure  5-8  shows  the  surface-focus  S-wave  SSSCs  for  station  CMTO  in  Thailand.  The 
correction  values  range  from  -  24. 1  seconds  to  +  41 .2  seconds.  The  values  are  near  zero 
at  the  station  and  show  greater  differences  as  we  get  to  greater  regional  distances.  The 
SSSCs  tend  to  have  zero  or  more  negative  values  for  locations  in  the  continental  areas  to 
the  north  and  from  the  shallow  oceanic  areas  directly  south.  Large  positive  SSSCs  occur 

for  locations  from  the  oceanic  areas  to  the  southwest  and  to  the  far  east.  These  again 
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appear  to  be  reasonable  for  these  regions  because  of  the  low  model  velocities  in  these 
oceanic  regions,  where  velocities  are  significantly  lower  than  the  IASPEI91  1-D  model 
velocities. 

5.4  DEPTH  DEPENDENCE  OF  SSSCs 

As  noted  above,  a  significant  advantage  of  our  model-based  approach  over  purely 
statistical  methods  for  predicting  SSSCs  is  that  we  can  develop  the  travel  time  corrections 
over  a  range  of  potential  focal  depths  for  sources  in  the  region.  For  each  of  the  30  IMS 
stations,  we  have  used  our  SSSC  prediction  scheme  to  predict  the  SSSCs  for  not  only 
focal  depths  at  the  earth’s  surface  but  also  for  10  additional  depths  within  the  crust  and 
upper  mantle  (0  km  <  h  <  200  km). 

Plots  of  the  SSSCs  at  stations  in  some  of  the  platform  regions  do  not  show  very  strong 
depth  dependence,  while  greater  dependencies  are  seen  at  several  of  the  stations  for 
events  with  sources  near  tectonic  boundaries  with  sharp  velocity  contrasts.  We  illustrate 
the  behavior  of  the  depth  dependence  of  the  SSSCs  for  sample  plots  at  selected  source 
depth  intervals  as  determined  for  the  P  waves  at  station  PRPK  at  Nilore  in  Pakistan.  The 
region  surrounding  station  PRPK  includes  some  fairly  strong  tectonic  boundaries,  such  as 
the  tectonic  plate  convergence  zone  in  the  Pamir-Hindu  Kush  area  just  to  the  north  of  the 
station. 

Figures  5-9  through  5-12  show  the  P-wave  SSSCs  at  PRPK  for  source  depths  of  0  km, 

10  km,  30  km,  and  100  km,  respectively.  The  P-wave  SSSCs  at  PRPK  for  regional 
surface  focus  sources  vary  from  -  7.3  seconds  to  +  5.3  seconds.  As  with  the  other  SSSC 
plots,  the  correction  values  for  source  locations  closest  to  the  station  are  small.  The 
correction  values  tend  to  be  negative  or  zero  over  large  parts  of  the  region  surrounding 
PRPK,  and  some  of  the  largest  negative  corrections  occur  out  near  the  regional  distance 
margins  (e.g.  strong  negative  anomaly  in  the  area  to  the  east  of  the  Caspian  Sea  in 
Turkmenistan).  Over  these  paths  the  3-D  model  travel  times  are  shorter  than  the 
corresponding  LASPEI91  times,  and  the  corresponding  velocities  in  our  3-D  model  are 
generally  higher  than  IASPEI91  1-D  model  velocities.  However,  a  remarkable  feature  of 
the  PRPK  SSSCs  in  Figure  5-9  is  the  strong  positive  corrections  (with  values  near  5 
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Figure  5-9.  Estimated  P-wave  SSSCs  at  the  PRPK  station  site  for  regional  events 
with  focal  depths  at  the  earth’s  surface  range  from  -7.3  seconds  to  +5.3  seconds. 
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events  with  focal  depths  of  10  km  range  from  -7.4  seconds  to  +5.1  seconds. 
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Figure  5-11.  Estimated  P-wave  SSSCs  at  the  PRPK  station  site  for  regional  events 
with  focal  depths  of  30  km  range  from  -7.2  seconds  to  +4.9  seconds. 
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Figure  5-12.  Estimated  P-wave  SSSCs  at  the  PRPK  station  site  for  regional  events 
with  focal  depths  of  100  km  range  from  -7.1  seconds  to  +3.6  seconds. 
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seconds)  to  the  north  in  the  Pamir-Hindu  Kush  region  of  Tadzhikistan.  For  these  regions 
the  travel  times  predicted  by  our  3-D  model  are  significantly  larger  than  the  IASPEI91 
times,  and  the  velocities  in  our  model  are  thus  much  slower  than  in  the  IASPEI91  model. 
Furthermore,  we  note  in  comparison  to  the  SSSC  plots  at  other  source  depths  that  there 
are  some  strong  variations  in  the  predicted  corrections  for  sources  in  the  Pamir-Hindu 
Kush  region  and  in  the  vicinity  of  the  India/Pakistan  coast.  Although  the  maximum  and 
minimum  SSSCs  are  about  the  same  in  the  plots  for  sources  with  focal  depths  of  10  km 
and  30  km  in  the  earth’s  crust,  there  are  some  significant  differences  in  the  correction 
values  at  specific  locations.  In  particular,  across  the  region  surrounding  PRPK 
differences  in  the  predicted  SSSCs  between  focal  depths  of  0  km  and  10  km  are  as  great 
as  2.5  seconds  in  some  places.  For  sources  with  subcrustal  depths,  the  positive  correction 
anomaly  in  the  Hindu  Kush  area  has  nearly  vanished.  In  general,  the  difference  in  the 
SSSCs  between  a  surface  focus  event  and  one  with  a  subcrustal  focal  depth  in  the  Hindu 
Kush  area  would  be  expected  to  be  as  great  as  2-3  seconds. 

Comparisons  like  those  in  Figures  5-9  through  5-12  indicate  that  application  of  SSSC 
correction  terms  from  a  fixed  or  incorrect  focal  depth  could  lead  to  substantial  errors  in 
the  travel  time  measurements  at  the  stations  where  they  are  applied.  Such  erroneous 
travel  times  can  lead  to  significant  location  errors  particularly  for  smaller  events  located 
with  limited  station  networks.  This  is  also  important  because  erroneous  travel  times  at 
stations  in  the  near-regional  distance  range  can  also  produce  significant  errors  in 
estimates  of  event  focal  depth  (cf.  Murphy  et  ah,  2000,  2002),  which  often  plays  a  key 
role  in  event  identification.  Therefore,  we  believe  that  additional  work  is  warranted  to 
establish  and  refine  SSSC  correction  and  location  procedures  to  account  for  these 
potential  effects  of  focal  depth  differences. 
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SECTION  6 


APPLICATION  OF  A  NEW  KRIGING  ALGORITHM  TO  ESTIMATE 
EMPIRICAL  SSSCs  AT  IMS  STATIONS 

6. 1  JOINT  LOCATION/CALIBRATION  INVERSE  PROBLEM 
Our  approach  to  seismic  calibration  determines  two  types  of  travel-time  functions:  a 
model-based  function  and  an  empirical  correction  function.  The  model-based  travel-time 
function  depends  on  three  variables:  the  event  location,  x,  the  station  location,  y,  and  a 
vector  of  model  parameters  that  describe  the  velocity  structure  of  the  Earth,  which  we 
will  denote  m.  Thus,  we  can  write  model-based  times  as  r(x,y;m).  Empirical 
corrections,  in  general,  depend  on  x,  y  and  a  parameter  vector  a:  C(x,y;a).  The  role  of  a 
is  analogous  to  that  of  m,  except  how  T  depends  on  m  is  determined  by  ray  theory  and  is 
not  a  matter  of  choice,  whereas  the  function  C  is  ad  hoc  and  so  the  definition  of  a  and 
how  C  depends  on  x,  y  and  a  are  a  matter  of  choice. 

The  problem  of  seismic  calibration  is  that  of  determining  m  and  a.  Our  approach  has 
been  to  find  both  vectors  by  solving  an  inverse  problem  involving  arrival  time  data  from 
multiple  events  and  stations,  including  where  possible  the  IMS  station  we  wish  to 
calibrate.  To  describe  our  approach,  we  consider  here  the  basic  version  of  this  problem 
with  first  P  arrival  times  observed  for  a  subset  of  paths  between  m  events  and  n  stations: 

dij  =  ti  +  T(xi, y/,m)  +  C(x,,y7;a)  +  e,y.  (6- 1 ) 

In  this  equation 

•  is  the  observed  arrival  time  for  the  (y')th  path. 

•  ti  and  X/  are  origin  parameters  (time  and  hypocenter,  respectively)  of  the  zth 
event. 

•  y j  is  the  location  of  the yth  station. 

•  e,y  is  an  observational  (picking)  error  in  dy. 

If  any  of  the  m  events  in  this  inverse  problem  are  earthquakes,  or  explosions  with 
imperfect  origin  information,  then  the  origin  parameters  of  these  events  are  additional 
unknowns  of  the  inverse  problem.  Therefore,  Equation  6-1  in  general  describes  a  joint 
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inverse  problem  whose  unknowns  are  m,  a,  and  (x„  /,•),  i  =  1, .. m.  As  part  of  this 
problem  we  include  appropriate  constraints  that  embody  whatever  ground-truth 
information  is  available  on  the  event  origin  parameters. 

We  have  solved  this  inverse  problem  geographic  area  by  area,  and  for  a  given  area,  in 
two  steps.  The  first  step  is  tomography,  in  which  a  is  fixed  (or  the  empirical  corrections 
are  ignored  altogether),  and  the  Earth  model,  m,  is  fit  to  the  data.  In  India-Pakistan,  we 
have  also  allowed  earthquake  location  parameters  to  vary  (i.e.,  earthquake  tomography). 

In  the  second  (post-tomography)  step,  we  solve  Equation  6-1  for  a  with  m  and  the  event 
locations  fixed  to  their  tomographic  results.  In  this  second  step,  it  is  useful  to  replace  the 
arrival  time  data  with  corresponding  post-tomography  residuals,  defined  as 

rg  =  dij  -  U  -  r(x„y7;m)  (6-2) 

and  write  the  inverse  problem  for  a  as 

ry  =  C(x„y7;a)+  etj.  (6-3) 

The  remainder  of  this  section  focuses  on  the  choice  of  a  and  the  estimation  of  a  from 
Equation  6-3. 

6.2  STATION-SPECIFIC  EMPIRICAL  CORRECTIONS 
6.2.1  PARAMETERIZATION 

The  traditional  approach  to  empirical  travel-time  corrections  defines  a  separate  correction 
function  for  each  station.  That  is,  there  is  no  functional  dependence  on  the  station 
location.  The  simplest  parameterization  of  station-specific  corrections  comprises  a  time 
term,  aj,  for  each  station,  such  that 


C(x,y7)  =  dj. 


(6-4) 


This  choice  eliminates  the  dependence  on  event  location  as  well,  which  is  not  an 
adequate  parameterization  for  spatially  well-distributed  events.  However,  it  is  the 
parameterization  used  in  many  multiple-event  location  methods  (e.g.,  Jordan  and 
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Sverdrup,  1981),  which  apply  to  small  clusters  of  events,  and  we  mention  it  for  this 
reason. 

Schultz  et  al.  (1998)  developed  an  empirical  calibration  approach  in  which  a  station- 
specific  correction  function  is  assigned  to  each  station: 


C(x,y;)  =  aj  (x). 


(6-5) 


In  this  case,  the  inverse  problem  of  equation  (6-3)  becomes 

ry  =  <2/(x/)+  ey  .  (6-6) 

This  problem  decouples  by  stations  since  no  connection  is  assumed  between  the 
correction  functions  for  different  stations.  For  any  given  station,  the  problem  is  one  of 
spatial  interpolation,  whereby  the  function  a/x)  is  fit  to  samples  observed  at  various 
points,  x,  (the  event  locations).  Schultz  et  al  (1998)  solved  this  interpolation  problem 
using  the  method  of  kriging.  Kriging  finds  an  optimal  fit  to  the  sample  values  under  prior 
constraints  in  the  form  of  statistical  properties  of  the  function  aj,  as  expressed  by  a 
correlation  function. 

6.2.2  APPLICATION  TO  POST-TOMOGRAPHY  DSS  RESIDUALS 
We  have  applied  station-specific  kriging  to  the  final  travel-time  residuals  resulting  from 
Pn  tomography  in  the  DSS  region  of  eastern  Asia.  Tomography  was  performed  in  two 
sub-areas:  one  centered  on  station  BRV,  which  we  will  refer  as  the  “western  DSS”  area, 
and  the  other  centered  35  degrees  longitude  east  of  BRV  (“eastern  DSS”).  We  applied 
kriging  to  the  post-tomography  residuals  from  each  area  separately,  even  though  many 
stations  were  involved  in  both  tomographies. 

We  used  the  2-D  “block  kriging”  program  KB2D  from  the  free  software  library  GSLIB 
(Deutsch  and  Joumel,  1998)  to  obtain  the  station-specific  travel-time  correction  function 
for  each  station.  Each  correction  function  is  2-D,  i.e.  a  function  of  geographic  position 
(event  epicenter)  but  not  event  depth.  Since  KB2D  works  in  a  Cartesian  geometry,  event 


latitudes  and  longitudes  were,  for  each  station,  converted  to  a  station-centered  Cartesian 
system.  The  KB2D  input  parameters  we  used  were 

♦  An  exponential  type  correlation  function:  cr02  exp(- |  x  |  Z) 

•  Correlation  distance:  L  -  400  km 

♦  Prior  standard  deviation:  =  1.5  sec 

•  Data  standard  deviations: 

-  0.5  sec  for  residuals  from  PNE  events 

-  1.0  sec  for  residuals  from  earthquakes 

We  show  the  station-specific  kriging  results  for  the  five  stations  for  which  the  most  data 
were  available:  BRV,  ELT,  NVS,  SVE  and  OBN.  Figure  6-1  shows  two  correction 
functions  for  station  BRV.  The  result  on  the  left  of  the  figure  used  the  227  residuals  at 
BRV  (from  48  PNE’s  and  179  earthquakes)  obtained  from  the  western  DSS  tomography. 
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Figure  6-1.  Correction  functions  for  station  BRV,  obtained  by  station-specific 
kriging  with  GSLIB.  A  correction  function  was  determined  from  the  residuals  from 
the  western  DSS  tomography  (left)  and  eastern  DSS  tomography  (right).  Locations 
of  the  events  in  each  data  set  are  shown  as  small  circles.  Five  stations  in  the  region 

are  shown  for  reference. 
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The  result  on  the  right  used  the  59  BRV  residuals  (46  PNEs,  13  earthquakes)  from  the 
eastern  DSS  tomography.  The  results  are  similar  except  the  western  DSS  correction 
function  (Figure  6-1  left)  has  an  area  of  large  negative  values  in  the  southeast,  not  present 
in  the  eastern  DSS  version  (right  of  figure).  This  is  probably  due  to  the  fact  that  the 
eastern  DSS  data  set  excluded  most  of  the  earthquakes,  some  of  which  had  negative 
residuals  as  large  as  -  5.0  sec. 

For  the  other  four  stations  we  show  only  the  correction  functions  obtained  from  the 
western  DSS  residuals.  Figure  6-2  shows  the  travel-time  correction  functions  for  ELT 
(left)  and  NVS  (right),  which  are  to  the  east  of  BRV.  The  correction  functions  for  these 
two  stations  are  similar,  which  is  expected  since  the  stations  are  separated  by  only  267 
km.  However,  it  is  difficult  to  know  whether  the  differences  that  we  do  see  (e.g. 
approximately  1  sec  near  station  SVE)  are  significant,  or  whether  they  are  artifacts  of  the 
different  event  sets  used  (35  PNEs  for  ELT,  26  PNEs  for  NVS)  or,  for  the  events  in 
common,  differences  between  the  observational  errors  at  the  two  stations.  Figure  6-3 
shows  the  kriged  correction  functions  for  stations  SVE  and  OBN,  two  stations  to  the  west 
of  BRV.  These  results  are  less  similar  to  each  other,  but  the  separation  between  SVE  and 
OBN  is  much  greater  than  ELT/NVS. 

Another  difficulty,  apparent  by  comparing  the  correction  functions  for  all  five  stations, 
has  to  do  with  source-receiver  reciprocity.  If  the  travel-time  corrections  embodied  in 
these  functions  are  the  result  of  Earth  structure  effects  absent  from  our  tomographic 
models,  the  correction  functions  should  be  source-receiver  reciprocal.  Thus,  the 
correction  predicted  for  station  ELT  from  a  shallow  event  located  where  NVS  is  should 
be  the  same  as  the  correction  for  NVS  from  an  event  located  at  ELT.  We  see  from  Figure 
6-2  that  this  is  nearly  the  case.  Comparing  Figures  6-1  (left)  and  6-2,  we  also  see  that 
BRV/NVS  reciprocity  is  obeyed  within  about  0.5  sec.  However,  BRV/ELT  reciprocity  is 
violated  by  more  than  1 .5  sec.  Reciprocity  between  BRV  and  SVE  is  also  violated  by 
more  that  1.5  sec  (Figure  6-1  versus  6-3).  Finally,  we  point  out  that  the  SVE  correction 
function  predicts  a  correction  of  about  -  1.5  sec  at  the  location  of  SVE  itself.  This  too  is 
non-physical. 


83 


Northing  (km) 


Kriged  Time  Residuals;  case  BRV  EA1  NEWc.P.ELT 
2000 


1500 
OBN 

ioooI-  ik 


g 

a> 


500 

0 


tr 

5  -500 


-1000 

-1500 

-2000 


^OSVE 

.  OBRV  °NVSt 


-2000  -1000  0  1000  2000 
Easting  (km) 


Kriged  Time  Residuals:  case  BRV  EA1  NEWc.P.NVS 
2000  ■ 


-2000  -1000  0  1000  2000 
Easting  (km) 


Figure  6-2.  Correction  functions  for  stations  ELT  (left)  and  NVS  (right),  each 
obtained  by  station-specific  kriging  of  the  western  DSS  residuals. 
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Kriged  Time  Residuals:  case  BRV  EA1  NEWc.P.OBN 


Figure  6-3.  Correction  functions  for  stations  SVE  (left)  and  OBN  (right),  each 
Obtained  by  station-specific  kriging  of  the  western  DSS  residuals. 
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Table  6-1  compares  statistics  of  the  travel-time  residuals  at  the  5  stations,  before  (prior) 
and  after  (posterior)  kriging.  We  see  that,  in  all  cases,  the  posterior  mean  residual  at  each 
station  (with  either  data  set)  is  essentially  zero,  and  the  posterior  standard  deviation 
(RMS)  has  been  reduced  significantly  from  its  prior  value. 


Table  6-1: 

Residual  Statistics  for  Station-Specific  Kriging  of  DSS  Residuals 

Station 

No.  Data 

Mean  (s) 

Std.  Deviation  (s') 

Prior 

Posterior 

Prior 

Posterior 

—Western  DSS— 

BRV 

227 

-0.059 

-0.023 

1.456 

0.629 

ELT 

35 

0.005 

0.003 

0.714 

0.326 

NVS 

26 

0.098 

-0.001 

1.070 

0.367 

SVE 

49 

-0.148 

0.003 

1.473 

0.494 

OBN 

29 

0.289 

-0.002 

0.840 

0.382 

— Eastern  DSS — 

BRV 

59 

0.138 

-0.038 

1.218 

0.513 

ELT 

35 

0.022 

0.003 

0.700 

0.330 

NVS 

26 

0.174 

0.003 

0.982 

0.320 

SVE 

11 

-0.179 

0.000 

1.597 

0.396 

6.3  MULTIPLE-STATION  EMPIRICAL  CORRECTIONS 
6.3.1  PARAMETERIZATION 

Rodi  et  al.  (2002, 2003)  have  developed  a  new  approach  to  empirical  calibration  that  fits 
correction  functions  to  a  multiple-event,  multiple-station  data  set,  such  as  is  used  in 
tomographic  model  development.  The  key  new  element  of  the  approach  is  the 
parameterization  of  path  corrections  with  a  universal  underlying  function  (or  set  of 
functions),  analogous  to  the  velocity  model  in  tomography.  Such  parameters  are  not 
specific  to  any  particular  station  or  event,  and  generate  a  correction  for  an  arbitrary  path. 


For  application  to  calibration  in  eastern  Asia,  we  have  considered  three  such 
parameterizations.  One  is  a  time-term  function,  a,  which  yields  the  correction  for  a  path 
between  event  location  x  and  station  locations  y  as 
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C(x,y)  =  a(x)  +  a(y). 


(6-7) 


The  parameter  vector,  a,  appearing  in  Equations  6-1  and  6-3,  would  in  this  case  be 
synonymous  with  the  function  a  or,  in  practice,  samples  of  this  function  on  a  dense  grid. 

It  is  clear  that  the  time-term  function  generates  reciprocal  travel-time  corrections,  since  x 
and  y  can  be  interchanged  in  Equation  6-7. 

Another  parameterization  we  have  considered  is  a  time-factor  function,  whereby 

C(x,y)  =  [a(x)  +  fl(y)]  T(x, y)  (6-8) 

where  T(x, y)  is  the  model-based  travel-time  for  the  path.  In  this  case,  the  correction  is 
proportional  to  the  travel-time  itself.  This  parameterization  also  yields  reciprocal  path 
corrections,  and  it  obeys  another  physical  property  of  travel-times  in  that  the  correction 
for  a  path  goes  to  zero  as  its  length  goes  to  zero: 

lim 

i  |  nC(x,y)  =  0.  (6-9) 

|x  -  y|  — >  0 

The  third,  and  last,  parameterization  we  have  considered  is  a  mislocation-vector  function, 
a(x),  for  which 


C(x,y)  =  p(x,y)  •  a(x)  +  q(x,y)  •  a(y).  (6-10) 

Here,  p  is  the  gradient  of  the  model-based  travel-time  with  respect  to  the  event  location, 
while  q  is  the  gradient  with  respect  to  the  station  location: 

p  =  V,r(x,y)  (6-11) 

q  =  V2T(x,y).  (6-12) 

(V4  means  gradient  with  respect  the  Ath  argument  of  a  function.)  Each  is  a  slowness 
vector  pointed  in  the  direction  of  the  raypath,  and  each  depends  on  both  x  and  y  since  the 
raypath  depends  on  the  location  of  both  its  endpoints.  The  interpretation  of  a  as  a 
mislocation  vector  follows  from  the  fact  that  Equation  6-10  is  the  first  order  change  in 
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travel-time  that  results  from  perturbing  the  event  location  by  a(x)  and  the  station  location 
by  a(y). 

Since  T  is  a  reciprocal  function  of  its  arguments,  we  have 

p(x,y)  =  q(y,x).  (6-13) 

Therefore,  the  mislocation-vector  parametrization  also  generates  travel-time  corrections 
that  are  source-receiver  reciprocal.  Further,  as  the  raypath  shrinks  to  zero  the  correction 
goes  to  zero  since  p  +  q  goes  to  zero. 

Station-specific  correction  functions  can  be  defined  from  these  universal  parameter 
functions.  In  the  case  of  Equation  6-10,  for  example,  they  are  simply 

aj(x)  =  p(x,y y)  •  a(x)  +  q(x,y7)  •  a(yy).  (6-14) 

This  points  up  an  important  advantage  of  universal  parameterizations  such  as  those 
above.  They  define  a  correction  function  for  a  station  at  any  location,  whether  or  not  the 
station  has  observed  travel-time  residuals  to  constrain  the  function.  This  is  accomplished 
by  allowing  correlations  between  the  corrections  among  stations  and  between  stations 
and  events. 

6.3.2  MAXIMUM-LIKELIHOOD  KRIGING 

The  estimation  of  universal  parameter  functions  from  observed  travel-time  residuals  is  an 
inverse  problem  which  is  a  generalization  of  the  interpolation  problem  encountered  in  the 
station-specific  approach.  For  example,  with  the  time-term  parameter  function,  the 
inverse  problem  of  Equation  6-3  becomes 

rij  =  a(x,)  +  a(yj)  +  ey.  (6-15) 

This  problem  (and  those  for  the  other  parameter  functions)  can  be  solved  with  geo- 
statistical  constraints,  as  used  in  conventional  kriging.  We  do  so  by  adopting  a 
maximum-likelihood  (ML)  approach  to  the  inverse  problem,  instead  of  the  minimum- 
variance  approach  used  in  conventional  kriging.  For  Gaussian  data  errors,  the  ML 
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approach  takes  the  solution  of  Equation  6-15  to  be  the  function  a  that  minimizes  the 
objective  function  given  by 

'i/  =  Zlrtf-fl(xj)-fl(yy)|2  /crJ+  Ja(x)[Dfl](x)cfx  (6-16) 

u 

where  erg  is  the  assigned  standard  deviation  of  eg  and  D  is  a  specified  differential 
operator.  This  is  very  similar  to  our  approach  to  solving  the  tomography  problem,  but 
now  we  choose  D  differently.  The  maximum-likelihood  inversion  approach  is  equivalent 
to  conventional  kriging  with  a  given  correlation  function  when  D  is  the  inverse  of  the 
covariance  operator  that  corresponds  to  the  correlation  function  (see  Rodi  et  al,  2003). 

In  our  applications  to  date,  we  have  used  the  following  operator  for  kriging  in  N 
spatial  dimensions  (N=  1,  2  or  3): 


D  = 


const 


<?(*)“ 


1 


2  l-N 


,2  &  ,2  &  ,2  & 

L  ~  +  Z  t  +  L  t 

.  y  dy1  1  dz1 J 


(6-17) 


where  o\  is  a  prior  variance  of  the  parameter  function  and  Lx,  Ly  and  Lz  are  correlation 
lengths  in  the  x,  y  and  z  directions  of  a  local  Cartesian  coordinate  system.  Various 
choices  of  the  operator  order,  i,  correspond  to  correlation  functions  that  are  commonly 
used  in  conventional  kriging,  e.g.  the  exponential  correlation  function  (£  =  1  in  1-D,  £ 
=  2  in  3-D)  and  Gaussian  correlation  function  (£  — >  oo).  While  we  have  expressed  this 
operator  in  Cartesian  coordinates,  we  have  actually  implemented  it  as  a  difference 
operator  in  a  global,  spherical  coordinate  system.  We  solve  the  maximum-likelihood 
criterion  (T  =  minimum)  for  a  sampled  version  of  the  function  a  on  a 
latitude/longitude/depth  grid.  Numerically,  the  minimization  is  performed  with  a 
conjugate  gradients  technique. 


6.3.3  APPLICATION  TO  POST-TOMOGRAPHY  DSS  RESIDUALS 
We  applied  the  multiple-station  kriging  approach  to  the  travel-time  residuals  resulting 
from  the  two  DSS  tomographic  applications.  This  time,  unlike  the  station-specific 
application,  we  pooled  the  data  from  the  two  sub-areas  (western  and  eastern  DSS), 
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yielding  a  data  set  with  1026  travel-time  residuals  from  261  events  and  70  stations.  We 
applied  the  multiple-station  kriging  method  to  these  data  to  estimate  parameter  functions 
of  the  three  types  described  above:  time-term,  time-factor  and  mislocation- vector 
functions.  We  estimated  2-D  (no  depth  dependence)  and  3-D  versions  of  each.  The  geo- 
statistical  parameters  we  used  were  similar  to  those  used  for  station-specific  kriging: 

•  Correlation  order:  £  =2  (exponential  correlation  in  3-D). 

•  Lateral  correlation  distance:  L  x=  Ly  =  400  km. 

•  Vertical  correlation  distance:  Lz  =  40  km  (for  3-D  functions). 

•  Prior  parameter  errors: 

-  Time  terms:  do =  1  -5  sec 

-  Time  factors:  ao=  0.01 

-  Mislocation  vectors:  cto=  5  km 

These  prior  errors  are  roughly  equivalent  when  mapped  to  travel-time  corrections.  The 
picking  errors  we  assumed  were  the  same  as  before:  0.5  sec  for  PNE  data  and  1 .0  sec  for 
earthquake  data. 

Figure  6-4  shows  the  2-D  time-term  (top)  and  time-factor  (bottom)  parameter  functions 
derived  with  multiple-station  kriging  of  the  DSS  residuals.  We  recall  that  a  travel-time 
correction  for  any  path  is  defined  from  the  time-term  function  (top  of  figure)  by  summing 
the  value  of  the  function  sampled  at  the  event  and  station  locations.  The  time-factor 
function  (bottom)  is  similarly  sampled  at  two  points,  but  the  sum  is  multiplied  by  the 
model-based  travel-time  to  obtain  a  travel-time  correction.  We  see  that  these  two 
parameter  functions  have  a  similar  spatial  dependence  when  fit  to  the  DSS  residuals. 
Figure  6-5  shows  the  three  components  of  the  2-D  mislocation-vector  function.  These 
define  a  travel-time  correction  for  a  path  in  terms  of  the  first-order  travel-time 
perturbation  that  results  by  moving  the  end-points  of  the  path  by  the  amount  and  direction 
indicated  by  the  mislocation  vector. 

Figures  6-6  and  6-7  through  6-9  show  3-D  versions,  at  three  depths,  of  the  time-term  and 
mislocation-vector  parameter  functions  derived  from  the  DSS  residuals. 
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Figure  6- 4.  2-D  time-term  (top)  and  time-factor  (bottom)  parameter  functions  for 
DSS  data  set  and  large  circles  mark  the  station  locations  (five  are  labeled). 
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Figure  6-5.  2-D  mislocation- vector  function  for  DSS:  north  (top),  east  (middle), 
and  depth  (bottom)  components. 


91 


TimeTerm  (sec)  Russia-BRV-3d1-a  Z-Okm 
30*  E  60*  E  90*  E  120*E  150*E 


Time  Term  (sec)  Russia-BRV-3d1-a  Z  -  50  km 


70*  N 
60*  N 


TimeTerm  (sec)  Russia-BRV-3d1-a  Z-  100km 


50*  N 
40*  N 
30*  N 


-2  -1.5  -1  -0-5  0 


Figure  6-6.  3-D  time-term  function  for  DSS,  at  depths  of  0  km  (top),  50  km  (middle), 
andTO'O'km  (Bottom). 
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Figure  6-7.  3-D  mislocation-vector  function  for  DSS:  north  component  at  depths  of 
0  km  (top),  50  km  (middle),  and  100  km  (bottom). 
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Figure  6-8.  3-D  mislocation-vector  function  for  DSS:  east  component  at  depths  of 
0  km  (top),  50  km  (middle),  and  100  km  (bottom). 
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Figure  6-9.  3-D  mislocation-vector  function  for  DSS:  depth  component  at  depths  of 
0  km  (top),  50  km  (middle),  and  100  km  (bottom). 
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Table  6-2  shows  the  posterior  standard  deviations  (RMS)  achieved  by  each  of  the 
parameter  functions.  These  can  be  compared  to  the  prior  RMS  for  the  pooled  DSS 
residuals,  which  is  1 .223  sec.  The  time-term  and  time-factor  functions,  both  the  2-D  and 
3-D  versions,  yield  relatively  large  posterior  RMSs,  especially  when  compared  to  station- 
specific  kriging  (Table  6-1).  This  is  not  unexpected  because,  in  these  cases  of  multiple- 
station  kriging,  there  is  one  parameter  function  fitting  the  data  from  all  stations,  rather 
than  one  per  station  as  in  station-specific  kriging.  Even  when  the  universal  function  is  3- 
D,  there  are  insufficient  degrees  of  freedom  to  reduce  variances  significantly.  The 
mislocation-vector  parameterization,  on  the  other  hand,  fits  three  functions  to  the  data 
(the  vector  components)  and  we  see  from  Table  6-2  that  this  provides  enough  degrees  of 
freedom  to  achieve  a  significant  RMS  reduction.  While  still  not  as  good  as  the  station- 
specific  case,  the  multiple-station  results  are  presumably  much  more  stable  by  having 
fewer  free  parameters. 

Table  6-2:  Posterior  RMSs  for  Multiple-Station  Kriging  of  DSS  Residuals 

Prior  RMS  =  1 .223  sec;  Number  of  data  =  1-26 

Time-Term  Time-Factor  Misloc.- Vector 

2- D  0.934  0.932  0.791 

3- D  0.904  0.898  0.692 


As  noted  earlier,  each  parameter  function  we  have  derived  defines  a  travel-time 
correction  function  for  each  station,  i.e.  the  correction  for  all  event  locations  and  a  fixed 
station  location  (see  Equation  6-14)  We  show  these  for  station  BRV,  as  derived  from 
four  of  the  parameterizations,  in  Figure  6-10  We  see  that  the  different  parameterizations, 
while  they  yield  varying  posterior  RMSs,  are  not  distinctly  different.  Apparently,  each 
captures  the  key  features  of  the  DSS  residual  patterns  for  BRV,  but  with  differing 
fidelity.  Figure  6-11  compares  the  2-D  mislocation-vector  version  to  the  station-specific 
function  for  BRV.  The  one  notable  difference  is  that  the  area  of  large  negative 
corrections  southeast  of  BRV  is  absent  in  the  multiple-station  function  (left).  This 
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Figure  6-10.  Travel-time  correction  functions  for  station  BRV,  generated  from  four 
different  parameter  functions:  2-D  time-term  function  (top  left),  2-D  time  factor 
function  (top  right),  2-D  mislocation- vector  function  (bottom  left),  and  3-D 
mislocation-vector  function  (bottom  right). 
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Figure  6-11.  Comparison  of  travel-time  correction  functions  at  BRV  obtained  from 
multiple-station  (left)  and  station-specific  (right)  kriging.  The  multiple-station 
function  is  based  on  the  2-D  mislocation-vector  parameter  function. 

suggests  that  this  part  of  the  single-station  function  for  BRV  was  either  not  well- 
determined  by  the  BRV  residuals  or  was  inconsistent  with  the  pattern  of  residuals  at  other 
stations. 

Finally,  Figure  6-12  shows  the  travel-time  correction  functions  for  the  other  four  stations 
considered  earlier.  Each  is  derived  from  the  2-D  mislocation-vector  parameter  function. 
Unlike  in  the  station-specific  versions,  stations  ELT  and  NVS  (top  panels  of  figure)  now 
have  nearly  identical  travel-time  corrections,  consistent  with  their  small  separation  (267 
km).  SVE  and  OBN  (bottom  panels)  are  less  correlated  since  they  have  a  greater 
separation.  By  design,  source-receiver  reciprocity  is  now  obeyed  between  all  station 
pairs. 
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Figure  6-12.  Travel-time  correction  functions  for  4  stations,  all  generated  from  2-D 
mislocation-vector  parameter  function:  station  ELT  (top  left),  station  NVS  (top 
right),  station  SVE  (bottom  left),  and  station  OBN  (bottom  right). 
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SECTION  7 


SUMMARY  AND  CONCLUSIONS 
7.1  PRINCIPAL  FINDINGS 

This  report  describes  the  results  of  a  three-year  project  designed  to  improve  capabilities 
to  locate  seismic  events  in  eastern  Asia  by  using  regional  calibration  of  seismic  travel 
times.  The  main  objective  of  this  research  project  has  been  to  calibrate  the  travel  time 
characteristics  of  the  earth’s  crust  and  upper  mantle  beneath  this  region  over  the  depth 
range  sampled  by  regional  seismic  ray  paths  to  each  of  30  IMS  stations  from  all  source 
locations  of  potential  interest.  The  research  reported  here  was  performed  by  a 
Consortium  of  institutions  led  by  SAIC  and  including  scientists  from  Massachusetts 
Institute  of  Technology’s  Earth  Research  Laboratory,  Weston  Geophysical  Corporation, 
and  the  Russian  Institute  for  Dynamics  of  the  Geospheres. 

One  of  the  principal  results  of  this  research  has  been  the  demonstration  that  a  3-D 
velocity  model  representative  of  this  Group  1  region,  based  on  prior  knowledge  and 
refined  using  GT  calibration  data,  can  be  used  to  determine  very  accurate  regional  travel 
times  for  P  and  S  waves.  Such  accurate  travel  time  predictions  enable  implementation  of 
source  specific  station  corrections  which  improve  the  accuracy  of  seismic  locations, 
particularly  for  smaller  events  which  may  be  located  with  observations  from  limited 
numbers  of  stations  and  often  include  data  from  complex  or  secondary  seismic  phases  at 
regional  stations.  The  travel  time  predictions  at  each  station  are  represented  as  a 
combination  of  model-based  travel  times  plus  empirical  corrections.  For  the  model-based 
component,  we  have  utilized  and  demonstrated  the  capabilities  of  a  state-of-the-art  3-D 
raytracing  program.  The  latter  provides  efficient  and  accurate  estimates  of  seismic  travel 
times  through  the  complex  velocity  structures  present  in  this  region.  This  raytracing 
algorithm  has  also  proven  to  be  valuable  in  computing  sensitivities  of  travel  times  to 
velocity  changes  throughout  the  model,  which  have  been  used  in  a  tomographic  inversion 
to  refine  and  improve  the  3-D  velocity  model  for  subregions  within  Group  1  which  have 
more  complete  calibration  data. 
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As  a  part  of  the  calibration  and  validation  efforts  for  the  Group  1  region,  a  large  database 
of  travel  time  observations  from  GT  events  has  been  assembled  and  archived.  These  data 
include  an  unique  FSU  explosion  database  including  GT  nuclear  tests  and  PNE 
explosions  covering  large  parts  of  the  region.  Travel  times  from  other  GT  events  for  the 
India/Pakistan  and  China/Mongolia  regions  have  also  been  included  in  this  calibration 
database  archive.  These  have  all  provided  important  data  for  refining  the  3-D  velocity 
models  in  each  of  the  sub-areas  which  make  up  the  overall  Group  1  region  of  eastern 
Asia. 

Validation  efforts  for  our  final  3-D  velocity  model  for  the  Group  1  region  have  included 
several  different  comparisons  and  analyses  of  seismic  travel  time  residuals  and  location 
errors.  These  validation  analyses  have  shown  that  P-wave  travel  time  errors  are  very 
small  for  our  refined  3-D  model  compared  to  nominal  IASPEI91  travel  times.  Therefore, 
the  corrected  travel  times  from  these  3-D  models  have  very  small  bias  and  the  RMS 
scatter  is  greatly  reduced.  Furthermore,  where  the  GT  events  have  been  available  to  test 
location  accuracy  in  the  Group  1  region,  we  have  found  significant  reductions  in  the 
location  errors  for  errors  for  events  located  with  the  3-D  model  travel  times  over  the 
IASPEI91  times. 

As  a  part  of  this  research  program,  we  have  also  conducted  preliminary  investigations 
using  this  same  approach  to  calibrate  travel  times  for  secondary  seismic  phases  from  the 
Group  1  region.  Although  this  approach  appears  promising,  some  additional  study  will 
be  needed  to  fully  demonstrate  its  potential  for  routine  event  location.  We  used  an 
analysis  of  Poisson’s  ratio  in  the  earth’s  crust  and  upper  mantle  to  determine  a  3-D  S- 
wave  velocity  model,  which  is  based  on  our  final  3-D  P-wave  velocity  model  for  the 
Group  1  region.  We  conducted  limited  testing  of  this  S-wave  model  and  the  associated 
travel  time  calculations  following  similar  raytracing  procedures  to  those  used  for  P. 

These  indicate  that  the  computed  travel  times  are  generally  consistent  with  S-wave  travel 
time  observations  from  GT  events  in  the  region.  However,  we  have  also  found  that  the 
travel  time  variance  for  the  corrected  S  waves  are  much  larger  than  for  the  corresponding 
P  waves.  This  suggests  that  for  even  a  very  good  S-wave  velocity  model,  the  S  times 
should  be  significantly  down-weighted  in  the  seismic  location  scheme. 
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Finally,  as  a  product  of  this  research  program,  the  final  P-  and  S-wave  SSSCs  for  the  30 
IMS  stations  in  the  Group  1  region  have  been  determined  and  delivered.  These  SSSCs 
have  been  provided  for  distances  out  to  20°  from  each  of  the  stations  and  for  a  range  of 
1 1  potential  focal  depths,  from  the  earth’s  surface  down  to  200  km.  Application  of  these 
SSSCs  to  regional  travel  time  observations  at  these  Group  1  stations  is  expected  to  result 
in  accurate  travel  times  which  can  be  used  with  observations  from  teleseismic  stations  to 
provide  better  locations  for  small  events  from  this  region  of  eastern  Asia. 

7.2  SOME  SPECIFIC  RESULTS 

In  addition  to  the  main  conclusions  and  deliverables,  this  Group  1  project  has  produced  a 
number  of  additional  results  which  have  contributed  to  improved  understanding  of 
seismic  location  calibration  across  this  region  of  interest  in  eastern  Asia.  Some  additional 
accomplishments  of  this  program  are: 

®  A  final  3-D  velocity  model  of  the  entire  Group  1  area  covering  the  region  of  eastern 
Asia  has  been  formulated  which  incorporates  both  P  and  S  velocity  distributions  at 
resolutions  ranging  from  1/3  to  1  degree  across  the  area  of  interest. 

•  A  version  of  the  Soviet  explosion  ground  truth  database  consisting  of  1000  carefully 
validated,  regional  distance  initial  P-wave  arrival  time  observations  from  distinct  FSU 
PNE  and  weapons  tests  has  been  completed  and  delivered. 

•  A  sophisticated,  fully  nonlinear  tomographic  inversion  module  has  been  formulated, 
implemented,  and  applied  to  the  refinement  of  the  upper  mantle  P-wave  velocity 
distributions  for  the  DSS,  China/Mongolia,  and  WINPAK3D  sub-regions,  where  the 
calibration  data  were  more  abundant.  These  sub-models  have  been  incorporated  into 
the  overall  3-D  model  for  the  larger  study  region. 

®  Numerous  features  found  in  the  tomographic  refinement  of  the  3-D  velocity  models 
have  been  seen  to  correlate  with  mapped  tectonic  structures  from  geologic  maps  of  the 
eastern  Asia  region.  This  supports  qualitatively  our  conclusion  that  the  tomographic 
inversion  procedure  is  producing  reasonable  results. 

®  Further  extensive  testing  of  the  tomographically  refined  P  velocity  models  using  the 
available  data  from  ground  truth  events  indicates  that  the  corresponding  travel  time 
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predictions  have  essentially  zero  bias  and  total  RMS  model  errors  on  the  order  of  1 
second  across  much  of  the  Group  1  study  region. 

•  The  estimated  errors  computed  for  the  refined  3-D  regional  velocity  model  represent 
significant  improvement  over  the  default  IASPEI91  model  in  the  FSU  region  which 
produces  an  average  bias  of  -  3.65  seconds  and  RMS  errors  of  more  than  2  seconds. 

•  Regional  seismic  event  location  accuracies  for  a  widely  distributed  sample  of  21  FSU 
PNEs  with  GTO,  using  observations  at  from  4  to  6  IMS  or  surrogate  stations,  show  an 
average  mislocation  of  only  6.9  km  using  the  refined  regional  model,  as  compared  to 
an  average  mislocation  of  20.4  km  for  the  same  data  using  the  LA.SPEI91  model. 

•  Similar  analyses  of  20  Novaya  Zemlya  explosions,  with  a  small  network  of  IMS  or 
surrogate  stations,  also  indicates  an  average  mislocation  less  than  10  km  for  the  3-D 
regional  velocity  model  compared  to  an  average  mislocation  of  44  km  obtained  for  the 
IASPEI91  model. 

•  Ground  truth  data  from  Chinese  National  Network  stations  for  Lop  Nor  explosions  and 
selected  better  located  events  from  the  India/Pakistan  region  have  produced  similar 
validation  results  supporting  our  final  refined  3-D  regional  velocity  models  for  the 
China/Mongolia  and  India/Pakistan  subregions. 

•  Knowledge  of  Poisson’s  ratio  in  the  earth’s  crust  and  upper  mantle  can  be  used  to 
convert  the  3-D  regional  P-wave  velocity  model  to  a  reasonable  S-wave  velocity 
model  which  can  be  used  to  predict  travel  times  of  secondary  seismic  phases. 

•  Measurement  errors  and  travel  time  prediction  uncertainties  for  regional  S  waves  and 
other  secondary  phases  appear  to  diminish  the  value  of  such  observations  for  event 
location. 

•  A  new  reciprocal  kriging  method  for  handling  empirical  travel  time  errors,  left  over 
after  modeling  corrections,  appears  promising.  Application  of  these  empirical 
corrections  on  top  of  the  model-based  corrections  results  in  reduced  RMS  errors  on  the 
order  of  0.7  seconds  in  regions  of  the  FSU  where  the  best  calibration  data  are  available 
for  testing  the  methods. 

•  The  small  residual  errors  which  are  left  after  the  combined  model-based  and  empirical 
corrections  appear  to  be  approaching  the  limits  of  improvements  to  be  expected  from 
seismic  travel  time  calibration  methods. 
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7.3  OPPORTUNITIES  FOR  FUTURE  DEVELOPMENT 

The  procedures,  velocity  models,  and  ground  truth  data  developed  under  this  project  are 
valuable  resources  for  use  in  locating  seismic  events  occurring  throughout  eastern  Asia. 
Application  of  the  SSSCs  developed  here  have  been  demonstrated  to  significantly  reduce 
location  errors  and  uncertainties  for  regional  seismic  events  determined  by  observations 
from  the  relevant  stations.  However,  there  are  a  number  of  areas  where  these  procedures 
and  results  could  be  improved. 

Additional  work  is  needed  to  improve  the  velocity  models  in  several  areas.  In  particular, 
information  on  the  regional  velocity  structure  from  several  of  the  subregions  (e.g.  far 
eastern  Russia,  southeastern  Asia,  and  adjacent  oceanic  regions)  included  in  this  large 
area  of  eastern  Asia  is  very  sparse.  Throughout  this  project  we  saw  evidence  that  the  use 
of  velocity  models  based  on  analyses  of  longer  period  surface  waves,  to  provide 
background  or  fill-in  where  more  complete  knowledge  was  lacking,  produced 
insufficiently  accurate  estimates  of  short-period  regional  body-wave  seismic  travel  times. 
More  complete  knowledge  based  on  regional  travel  times  would  be  most  useful  to 
improve  the  models  in  several  of  these  subregions. 

Equally  important  to  better  knowledge  of  velocity  models  is  expanded  calibration  data. 
The  most  reliable  travel  time  predictions  can  be  generated  in  regions  where  the  velocity 
models  have  been  tested  and  validated  with  good  observations  from  accurate  ground  truth 
events.  We  found  in  our  tomographic  analyses  that  the  velocity  models  could  be  refined 
and  improved  in  several  parts  of  the  eastern  Asia  study  region  (e.g.  DSS,  WINPAK,  and 
China/Mongolia)  where  calibration  data  were  more  complete.  Additional  ground  truth 
events  and  calibration  data  could  provide  better  velocity  models  in  many  of  the 
subregions  of  eastern  Asia.  Such  data  are  particularly  needed  in  some  of  the  areas  (such 
as  those  cited  above)  where  there  is  little  prior  knowledge  available  of  the  velocity 
structure.  Even  some  of  the  areas  where  tomographic  analyses  were  conducted  during 
this  project  could  be  further  refined  using  additional  and  more  reliable  calibration  data  to 
fill  in  some  coverage  gaps.  For  example,  more  accurate  ground  truth  from  selected 
events  in  the  China/Mongolia  region  could  lead  to  additional  tomographic  refinements 
and  model  improvements  over  parts  of  that  subregion. 
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During  the  course  of  this  project,  work  on  travel  time  calibration  for  secondary  seismic 
phases  was  very  limited.  Although  the  eastern  Asia  model  which  we  developed  for  S- 
wave  velocities  looks  very  promising,  we  performed  only  very  limited  validation  tests. 
Additional  ground  truth  and  accurate  S-wave  travel  time  observations  would  be  very 
valuable  for  improving  the  S  velocity  model  in  this  region.  As  noted  above,  uncertainties 
in  travel  time  picks  for  S  waves  and  other  secondary  phases  are  particularly  problematic 
for  utilizing  these  signals  in  location  schemes.  Thus,  an  important  element  for  improving 
utilization  of  these  secondary  seismic  phases  in  regional  location  will  be  development  of 
consistent  timing  procedures.  This  is  especially  a  problem  for  the  explosion  sources,  for 
which  S  signal  onsets  are  often  quite  emergent.  Unfortunately,  it  is  just  such  explosion 
sources  which  most  often  provide  the  best  ground  truth.  So,  there  is  a  clear  need  for 
significant  additional  work  in  data  collection  and  analysis  if  secondary  seismic  phases  are 
to  be  made  useful  for  improving  the  location  accuracy  for  regional  seismic  events. 

Finally,  the  advanced  kriging  procedures  which  have  been  developed  as  part  of  this 
project  appear  to  provide  a  more  realistic  approach  to  represent  the  empirical  errors  and 
travel  time  uncertainties.  In  this  project  we  have  mainly  tested  these  procedures  in 
subregions  where  the  calibration  data  were  more  abundant.  More  complete  analyses  of 
these  reciprocal  kriging  methods  and  their  sensitivities  to  model  and  travel  time 
measurement  uncertainties  throughout  this  Group  1  study  region  as  well  as  other  regions 
of  interest  would  be  valuable. 


SECTION  8 


APPENDIX.  SURFACE  FOCUS  P-WAVE  SSSCs  AT  GROUP  1  IMS  STATIONS 

P-wave  and  S-wave  SSSCs  have  been  computed  using  the  procedures  based  on  three- 
dimensional  seismic  velocity  models  described  in  this  report  for  all  30  IMS  station  sites 
in  the  Group  1  calibration  region  of  eastern  Asia.  The  figures  on  the  following  pages 
show  plots  of  the  surface-focus  P-wave  SSSCs  which  illustrate  the  results  accomplished 
on  this  project.  The  digital  versions  of  these  SSSCs  along  with  those  from  10  additional 
focal  depths,  covering  potential  focal  depths  throughout  the  earth’s  crust  and  upper 
mantle  (h  <  200  km)  have  been  computed  and  delivered. 
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Figure  8-1 .  Surface  focus  P-wave  SSSCs  for  station  AAK  range  from 
-4.6  seconds  to  +5.4  seconds. 
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Figure  8-2.  Surface  focus  P-wave  SSSCs  for  station  AKTO  range  from 
-8.1  seconds  to  +3.1  seconds. 
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Figure  8-3.  Surface  focus  P-wave  SSSCs  for  station  BIL  range  from 
-3.8  seconds  to  +5.0  seconds. 
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Figure  8-4.  Surface  focus  P-wave  SSSCs  for  station  BJT  range  from 
-5.0  seconds  to  +4.3  seconds. 
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Figure  8-5.  Surface  focus  P-wave  SSSCs  for  station  BRVK  range 
from  -9.2  seconds  to  +3.0  seconds. 
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Figure  8-6.  Surface  focus  P-wave  SSSCs  for  station  CMTO  range  from 
-2.9  seconds  to  +5.8  seconds. 
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Figure  8-7.  Surface  focus  P-wave  SSSCs  for  station  EVN  range  from 
-4.0  seconds  to  +6.7  seconds. 


-4.4  seconds  to  +3.4  seconds. 
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Figure  8-9.  Surface  focus  P-wave  SSSCs  for  station  JUKI  range  from 
-3.1  seconds  to  +5.1  seconds. 
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Figure  8-10.  Surface  focus  P-wave  SSSCs  for  station  KMI  range  from 
-2.7  seconds  to  +4.9  seconds. 
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Figure  8-11.  Surface  focus  P-wave  SSSCs  for  station  KSRS  range  from 
-3.8  seconds  to  +7.0  seconds. 


Figure  8-12.  Surface  focus  P-wave  SSSCs  for  station  KURK  range  from 
-8.1  seconds  to  +3.8  seconds. 
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Figure  8-13.  Surface  focus  P-wave  SSSCs  for  station  LZDM  range  from 
-3.2  seconds  to  +3.3  seconds. 
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Figure  8-14.  Surface  focus  P-wave  SSSCs  for  station  MA2  range  from 
-3.8  seconds  to  +5.2  seconds. 
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Figure  8-15.  Surface  focus  P-wave  SSSCs  for  station  MKAR  range  from 
-6.4  seconds  to  +4.0  seconds. 
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Figure  8-16.  Surface  focus  P-wave  SSSCs  for  station  NRIK  range  from 

-8.9  seconds  to  +3.3  seconds. 
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Figure  8-17.  Surface  focus  P-wave  SSSCs  for  station  PALK  range  from 
-2.8  seconds  to  +6.5  seconds. 
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Figure  8-18.  Surface  focus  P-wave  SSSCs  for  station  PDY  range  from 
-8.8  seconds  to  +3.6  seconds. 
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Figure  8-19.  Surface  focus  P-wave  SSSCs  for  station  PRPK  range  from 
-7.3  seconds  to  +5.3  seconds. 


Figure  8-20.  Surface  focus  P-wave  SSSCs  for  station  SF.Y1  rangf^fmm 

-2.8  seconds  to  +6.7  seconds. 
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Figure  8-21 .  Surface  focus  P-wave  SSSCs  for  station  SONM  range  from 
-5.3  seconds  to  +2.7  seconds. 
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Figure  8-22.  Surface  focus  P-wave  SSSCs  for  station  SSE  range  from 
-3.7  seconds  to  +5.6  seconds. 
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Figure  8-23.  Surface  focus  P-wave  SSSCs  for  station  TIXI  range  from 
-8.0  seconds  to  +5.1  seconds. 
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Figure  8-24.  Surface  focus  P-wave  SSSCs  for  station  TLY  range  from 
-6.9  seconds  to  +3.0  seconds. 
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Figure  8-25.  Surface  focus  P-wave  SSSCs  for  station  URG  range  from 
-3.6  seconds  to  +6.0  seconds. 
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Figure  8-26.  Surface  focus  P-wave  SSSCs  for  station  USK  range  from 
-4.6  seconds  to  +7.3  seconds. 
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Figure  8-27.  Surface  focus  P-wave  SSSCs  for  station  XAN  range  from 
-4.1  seconds  to  +4.0  seconds. 
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-8.2  seconds  to  +5.4  seconds. 
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Figure  8-29.  Surface  focus  P-wave  SSSCs  for  station  YSS1  range  from 
-2.4  seconds  to  +7.3  seconds. 
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Figure  8-30.  Surface  focus  P-wave  SSSCs  for  station  ZAL  range  from 
-7.8  seconds  to  +3.4  seconds. 
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2  CYS  ATTN:  TDND/D.  BARBER 


DEPARTMENT  OF  DEFENSE  CONTRACTORS 

ITT  INDUSTRIES 
ITT  SYSTEMS  CORPORATION 
1680  TEXAS  STREET,  SE 
KIRTLAND  AFB,  NM  87117-5669 
2  CYS  ATTN:  DTRIAC 
ATTN:  DARE 

SCIENCE  APPLICATIONS  INTERNATIONAL 

CORPORATION 

10260  CAMPUS  DRIVE 

MAIL  STOP  XI 

SAN  DIEGO,  CA  92121 

ATTN:  J.R.  MURPHY 


